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ABSTRACT 



6 ? *3 £> 6 * 6 ? 

The Regenerator, Cooler and Heater for the NASA Space Power Research Engine ^ .i*7 


(SPRE) have been analyzed in detail for laminar, incompressible and oscillatory flow 
conditions. Each component has been analyzed independently and in detail with the 
regenerator being modeled as a^parall el -plates channel with a solid wall. The ends of the 

p Yo \A cb ^ 

channel are exposed to two reservoir maintained at different temperature thus fasdiitatihg 

(Am, 

^ axial temperature gradient along the channel. The cooler and heater components have 
been modeled as circular pipes with isothermal walls. Two different types of thermal 
boundary conditions have been investigated for the cooler and heater, namely, symmetric 
and asymmetric temperature inflow. In symmetric temperature inflow the flow enters the 
channel with the same temperature in throughout the velocity cycle whereas, in 
asymmetric temperature inflow the flow enters with a different temperature in each half 
cycle. The study was conducted over a wide range of Maximum Reynolds number 
(i?e max ) varying from 75 to 60000, Valensi number (Va) from 2.5 to 800, and relative 
amplitude of fluid displacement (A r ) from 0.357 to 1.34. 

A two dimensional Finite Volume method based on the SIMPLE algorithm was 
used to solve the governing partial differential equations. Post processing programs were 



developed to effectively describe the heat transfer mechanism under oscillatory flows. 
The computer code was validated by comparing with existing analytical solutions for 
oscillating flows. 

The thermal field have been studied with the help of temperature contour and three 
dimensional plots. The instantaneous friction factor, wall heat flux and heat transfer 
coefficient have been examined. It has been concluded that in general, the frictional 
factor and heat transfer coefficient are higher under oscillatory flow conditions when the 
Valensi number is high. Also, the thermal efficiency decreases for low'er A r values. 
Further, the usual steady state definition for the heat transfer coefficient does not seem 
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NOMENCLATURE 


Half channel width in the conjugate problem. 
Relative amplitude of fluid motion. 

Half solid thickness in the conjugate problem. 
Specific Heat of the fluid or solid. 

Hydraulic diameter of the tube or channel. 

Thermal conductivity. 

Length of the channel or tube. 

Nusselts’s number. 

Hydrodynamic pressure in the momentum equations. 
Peclet number used in the analytical solution. 
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Prandtl number of the fluid. 

Heat flux. 

Radial coordinate distance and direction for Axisymmetric cases or, 
the radius of the tube. 

Reynolds number based on the maximum velocity [«J, 

Time [s]. 

Temperature in °K,T(j:,y,0 or T{x,r,t) . 

Section averaged temperature of the fluid, T a {x,t ) 

JYrdr 

[rdr 

Bulk temperature of the fluid weighted by the absolute velocity, 

T b (x,t) 

Time averaged temperature, T x (x,y ) , 

+ . 

Axial velocity in the x coordinate direction, u(x,y,t) or u(x,r,t ) . 
Representative axial velocity for the analytical solution. 

Axial velocity at the inlet, u in (t) . 

Maximum input velocity in a cycle. 

Control volume. 

Normal velocity of the fluid in the y coordinate direction. 

Valensi Number. 


/ \u(r)\Trdr 
f\u{r)\rdr 



x Axial coordinate distance and direction. 

v Normal coordinate distance and direction for Cartesian Coordinate, 

system. 

SUBSCRIPTS 


c 

h 

i 

in 

in 

max 

fd 

w 


Cold end or east side. 

Based on hydraulic diameter or hot end reservoir. 
Instantaneous values. 

Inlet condition. 

Mean or averaged over the cross-section value. 

maximum value during the cycle. 

based on fully developed flow conditions. 

At the wall. 


SUPERSCRIPTS 


* Guessed or lagged quantity. 

GREEK SYMBOLS 


U) 

P 

Y 

v 

P 

o 

K 

e 

X 

a 

*1 


Angular frequency in rad/sec. 

Dynamic viscosity of the fluid. 

= -(dT/dx) = [T }vtst -T easr \ / L ,time averaged constant axial temperature- 
gradient. 

Kinematic viscosity of the fluid. 

Density of the fluid. 

Ratio of fluid to solid thermal diffusivity, 

=(k/pC p ) f Kkl P C p ) s . 

=kj k ,ratio of fluid to solid thermal conductivity. 

/ s 

= bt a , ratio of half solid to half fluid thickness. 

Non-dimensional Pressure gradient in the analytical solution, 

= \dp/dx\a 2 /pu t v. 

Womerseley number used in the analytical solution, 

=a\[w]\ . 

=y / a, nondimensional normal direction. 



CHAPTER I 


INTRODUCTION 

The Stirling engine is a efficient power producing device based on the Stirling 
thermodynamic cycle. A few of its distinct features include high efficiency, very long 
life, high reliability, and low noise. But, the most important feature which makes it a 
strong candidate for a viable power source in the future is its ability to driven by virtually 
any power source such as solar energy. Thus it is seen as a ideal power source for 
remote applications such as space systems and remote terrestrial applications. 

The Stirling Technology branch of NASA Lewis Research Center, Cleveland, have 
been working on developing free-piston Stirling engines for both space and civil 
applications (high and low power technology). Currently, research is going on to better 
understand the thermodynamic losses in the NASA SPRE (Space Power Research 
Engine), such that the efficiency can be maximized. A team from Cleveland State 
University have been working on two dimensional modelling and analysis of the SPRE 


z 



SPRE cran mc*w 


Figure 1.1 Cross sectional view of the NASA Stirling Space Power Research Engine. 
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heat exchangers to get a feel for the flow and heat transfer phenomenon occurring in the 
heat exchangers. The thermal analysis of the heat exchangers done in this study is based 
on the operating conditions of SPRE. Figure 1.1 shows a quarter sectional view of the 

NASA’s SPRE and its various components. 

Figure 1.2 is a schematic representation of the major components of the free-piston 
Stirling showing the basic components and their relative locations. The basic power 
output of the engine is based on the net work done on the piston by the working fluid or 
gas. The energy inputs include the heat input to the heater for heating the gas. This heat 
input is reduced by the addition of the regenerator thus making it a highly efficient 
thermodynamic cycle. The oscillatory motion is achieved by the pressure changes in the 
piston and displacer gas-springs thereby shuttling the working fluid to and fro from the 
compression to the expansion space. In the shuttling process the gas absorbs heat energy 
from the heater part of which is absorbed by the regenerator when the gas is on its way 
to the cooler. And when the gas flows from the cooler to the heater the regenerator 
releases this stored energy thereby reducing the net heat input to the cycle or engine. 

The location of the heat exchangers between the compression and expansion space 
results in an oscillatory flow in the heat exchangers. Therefore the design of the heat 
exchangers (heater, cooler or regenerator) needs to consider the effect of oscillatory flows 
on the thermal losses. Due to lack of sufficient data and analysis the heat exchangers 
designers still use unidirectional, steady state correlations for the friction factor and heat 
transfer coefficient. Since the flow losses and thermal losses work against each other the 
phenomenon has to be understood to optimize the engine working conditions. 
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Some analytical solutions have been derived for the flow between constant area 
ducts assuming fully developed flows. But the thermal analysis and analytical solution 
has been restricted to the presence of a linear axial temperature gradient (Kurzweg, 
1985a). No thorough analysis has been done on different temperature boundary 
conditions as present in the Stirling engine heat exchangers. 

1.1 Objectives of the research 


The present study concerns itself with the time dependent flow and thermal fields 
in the heat exchangers of the NASA SPRE, namely the, cooler, heater and regenerator. 
The flow in all the components was assumed to be laminar and incompressible. For the 
analysis effort was made to model the appropriate geometry for each of the heat 
exchangers and choose an efficient and reliable numerical method to solve for the 
governing equations needed for the analysis. Once the flow and thermal fields were 
established, the study focussed on the behavior of the instantaneous friction factor and 
heat transfer coefficient with the: 

► Maximum Flow Reynolds Number 

► Valensi Number or dimensionless frequency ( va) 

► Relative amplitude of fluid displacement (A r )- 

The study was conducted over a wide range of the above mentioned parameters 
in order to correctly assess the effect of these parameters on the flow and thermal field. 
Also, different thermal boundary conditions namely, the conjugate heat transfer type, 
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symmetric temperature inflow and asymmetric temperature inflow were also explored 
under oscillating flow conditions. Special efforts were made to validate the numerical 
method by comparing the predictions with existing analytical solutions. 
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CHAPTER II 

LITERATURE REVIEW 


The aim of this chapter is to synopsize the past work done on oscillatory flows. 
The survey discusses any analytical solutions and experimental data available on 
oscillatory flows. Also numerical simulations that have been done on oscillatory flows 
have been addressed and finally there is a discussion on further investigations that are 

needed to understand the effects of oscillatory flows. 

There are two kinds of unsteady (cyclic) flows that one can find in the literature, 
pulsatile and oscillatory flows. In pulsatile flows the fluid is set to motion by a 
sinusoidally varying pressure gradient or velocity which has a non zero mean, which 
means in a complete cycle there exists a net mass transfer across any cross section normal 
to the primary flow direction. The non zero mean also implies that the primary direction 
of the inflow does not change in a cycle. Whereas in oscillatory flows the flow is driven 
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by a harmonically varying pressure gradient or velocity that has a zero mean. Given the 
definitions of these two types of unsteady flows one is easily led to conclude that 
oscillatory flows are a special case of pulsatile flows when the driving pressure gradient 
or velocity approaches zero mean, which is not true at least mathematically since zero 
mean flow is a singularity point. This singularity makes the oscillatory flows phenomena 
a very complex one, although qualitatively one can observe similar behavior between the 
two kinds od flows. Since the present study concerns itself with only oscillatory flows 
the review of pulsatile flows has been left out and interested readers can find an extensive 
review of pulsatile flows in Kohler (1990) and Kw r an (1992). 

Before the survey is presented a brief description about the meaning of the term 
"oscillatory flows" in the present context needs to be elaborated. One can find in the 
literature about oscillatory flows in external flows (see Schlichting) where the flow pattern 
around a harmonically oscillating body immersed in a fluid are discussed. That situation 
is different from the one encountered in this study which is mainly concerned with 
internal flows. The characteristic of oscillatory flows in a internal flow situation is that 
the periodic driving force has a zero mean for a complete cycle, physically this means in 
a whole cycle there is no net mass transfer across any cross section perpendicular to the 
direction of the periodic input. Furthermore, in oscillatory flows because of the zero 
mean the direction of the flow is actually reversed from one half cycle to the other half 
cycle. 

As the effects of oscillatory flows are completely different from unidirectional 
flows which means the transition from laminar to turbulent oscillating flow is different. 
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Many experimental studies have been done on the transition criterion from laminar to 
turbulent flow in oscillatory flows and a comprehensive review can be found in Seume 
and Simon (1986a). Only a brief review of the experimental work directly related to the 
present study will be presented. 

2.1 Analytical Solutions 

One of the first analytical solution for the oscillatory flow problem was derived 
by Stokes, who obtained the flow field about an infinite flat wall which executes a 
sinusoidal motion ( See Schlichting ) in a stagnant fluid. The effect of the unsteady 
motion on the flow field was recognized by the presence of what is now known as the 
Stoke s layer. This layer is a small region close to the wall where the viscous diffusion 
is concentrated and the region away from it is not effected at all by the motion of the 
plate. Kurzweg and Chen (1988) did a heat transfer analysis on the above harmonically 
oscillating plate when it is subjected to a constant axial temperature gradient. 

Richardson (1928) in an acoustic experiment measured the velocity distributions 
across an orifice of circular cross section and he found the peak velocity close to the wall 
instead of the centerline of the orifice. This was theoretically verified by Sexl(1930) and 
experimentally corroborated by Richardson and Tyler (1929-30) for the flow produced by 
the reciprocating motion of a piston. They had mistakenly characterized the velocity 
peaks as "annular effect" due to the circular geometry. But it has been shown later these 
velocity shoots near the wall are characteristic of oscillatory flows even in parallel plates 
situation and not due to any particular geometry. 

From the literature surveys it appears that the laminar fully developed 
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oscillator}' flows are fairly well documented. The fully developed solutions were derived 
analytically by using parallel flow assumption and neglecting the initial conditions such 
that the flow develops into periodically steady state. All the analytical work has been 
done for constant area ducts and these are reviewed next including relevant numerical 
simulations and experimental studies. 

2.1.1 Two Dimensional Geometries 

Uchida (1956) calculated the velocity profiles for laminar incompressible flow in 
a circular tube subjected to a arbitrarily varying time dependent pressure gradient. He 
linearized the Navier-Stokes equation by assuming parallel flow thereby dropping the 
axial diffusion terms and was able to exactly solve the momentum equations by Fourier 
decomposition of the time dependent pressure gradient term. His paper also appears to 
be one of the first to distill out practically useful quantities as the wall shear stress. 

Kurzweg (1985a, 1985b) was one of the first to extend the analysis to include heat 
transfer for both parallel plate and circular geometry. It is appropriate to mention that all 
the analytical heat transfer solutions were derived for a thick walled 2D geometry, that 
is for the conjugated heat transfer problem. Kurzweg’s analysis was based on the earlier 
works of Chatwin (1975), Watson(1983), Joshi et al. (1983) whom found the diffusion 
of contaminants in gases were greatly enhanced when subjected to flow oscillations. 
Drawing an analogy for the diffusion of heat Kurzweg was able to arrive at a closed-form 
solutions for the temperature distribution in the channel. His findings indicate that in 
oscillatory flows if the fluid entering the channel has different specific enthalpy in one- 
half period than the other half, then the axial heat transfer is greatly enhanced due to fluid 
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oscillations. This heat transfer is further augmented if the channel wall has a finite 
thickness as it increases the heat storage/release capability of the wall. Kaviany (1990) 
and Zhao Ling-de et al. (1991) did a similar analyses for circular tubes with similar 
qualitative findings and in addition they also verified their theoretical predictions with 
experiments. It should be mentioned here that this enhanced axial heat diffusion is an 
undesirable effect for Stirling Engine application as the function of the regenerator (which 
is sandwiched between the Heater and Cooler maintained at different temperatures) is to 
minimize the axial transfer of heat from the Heater to Cooler and vice-versa. 

Gideon (1986) using a mean parameter approach (average over cross section) 
technique was able to arrive at practically useful relationships such as friction factor and 
heat transfer coefficient for the oscillating flow in a channel, but in general the results 
were not in agreement with the exact solutions. Although the technique proved useful for 
one dimensional modeling of the flow field in the Stirling Engine Heat Exchangers. 

Ozawa and Kawamoto (1991) correlated this enhanced axial diffusion in terms of 
an effective Nusselt Number for a range of Prandtl numbers ( p r ) and Reynolds Number. 
They observed that for higher p r the lateral diffusion of heat and momentum penetrated 
the same distance in other words the fluid behaved as though its Pr was unity for 
sinusoidal motion of fluid in a pipe. Based on this observation they used a lumped- 
parameter approach (two layered model) and arrived at a closed-form solution for the heat 
transfer coefficient also they validated these correlations experimentally. 

Kaviany (1986) extended the analytically investigations to include the effects of: 
(i) viscous dissipation, (ii) channel spacing ( height of the channel ), and (iii) wall 
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thickness (Solid Thickness). 

The above analytical investigations for laminar oscillating flows in constant area 
two dimensional ducts with parallel flow assumptions can be summarized as follows: 

- Additional external work is needed to sustain this type of flow due to 
higher friction losses as the wall shear stress is augmented. 

- When the frequency of oscillation is low the fully developed velocity 
profile approaches the parabolic shape found in steady, whereas at high 
frequencies the lateral momentum diffusion is restricted to the Stoke s 
layer. 

- Large quantities of heat can be transported across reservoirs, maintained 
at different temperatures connected by a channel without convective 
transport of mass (Oscillatory flows). 

- There exists a particular frequency when the axial heat transport is 
maximized. 

2.2 Numerical Simulations 

Ibrahim et al. (1989,1990) carried out numerical simulations assuming hydrody- 
namically developed flow and confirmed the analytical findings. Further more the effect 
of constant temperature boundary conditions on the heat transfer coefficient was also 
presented. Devalba M. et al. (1991) also simulated the conjugate heat transfer problem 
under oscillating flow conditions using a finite element code and confirmed the analytical 
predictions. 

Ahn (1990) and Kohler (1990) conducted extensive numerical investigations under 
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turbulent flow conditions utilizing various turbulence models. Their findings indicate that 
the two equation turbulence models are inadequate to correctly capture the effects of 
accelerating and decelerating flows. Further turbulent and transition studies in oscillating 
flows are documented in Koehler’s doctoral thesis. 

Recently Patankar and Oseid (1992) carried out a two dimensional heat transfer 
analysis in a pipe under turbulent flow conditions. The boundary conditions were similar 
to tha! of the Heater in the Stirling Engine and their study indicates the sensitivity of the 
Nusselt number (lateral) to the phase shift between the wall heat flux and the bulk 
temperature. 

Hashim (1992) studied the effect of fluid oscillations on the heat transfer and skin 
friction coefficient for various configurations of backward facing steps. Kwan (1992) 
investigated the compressibilty effects in a channel for oscillating and pulsatile flows. 

2.3 Experimental Studies 

A large part of the experimental research has been devoted to study the stability 
and transition mechanisms of oscillatory flows and can be found in Seume (1988). 
Seume (1988) carried out a number of experimental runs for different operating 
parameters in order to understand the transition mechanisms in oscillatory flows and 
eventually came out with an envelope identifying the laminar and fully turbulent regimes 
on a Re (Reynolds number based on the maximum velocity input and hydraulic 

111ILX 

diameter) and Va (Valensi number) plot. Friedman (1991) later made detailed 
3measurements at a particular operating point and with that database was able to extract 
useful information about the effect of fluid oscillations on wall bounded flows. 
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Yuan and Dybbs (1992) experimentally studied and simulated the regenerator of 
a Stirling engine under oscillatory flow conditions. The study was mainly concentrated 
on the effect of high frequency and high pressure on the heat transfer coefficient between 
the working fluid and solid matrix in the regenerator. Their findings indicate oscillation 
frequency effects both the temperature and heat transfer coefficient while the pressure 
effects only the heat transfer coefficient. Furthermore the heat transfer coefficient is 
enhanced significantly compared to that in unidirectional flows. 

2.4 Summary 

The survey presented above can be summarized as follows: 

- Laminar fully developed flows seems to be well understood for 
oscillatory flows. But the heat transfer analysis has been concentrated 
mainly on the constant heat flux (boundary condition) problem. 

- Currently there seems to no general consensus on the non dimensional 
parameters to be used especially when it comes to the dimensionless 
frequency. It has been referred to as Kinetic Reynolds number, Valensi 
Number and Womersely number. But the trend among the Stirling engine 
researchers seems to be adopt the definition of Valensi number. 

- Because of the fluid oscillations the velocity profile and hence the 
temperature profile take different shape as compared to unidirectional 
flow except at low oscillation frequencies. 

- In genera] all the physical quantities such as the pres- 
sure, velocities, temperature etc. are out of phase relative to each 
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significantly alter the friction factor and heat transfer coefficient. 

- As far as laminar to turbulence transition is concerned the respective 
regions are fairly well charted out and documented. 
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CHAPTER m 

MATHEMATICAL DESCRIPTION OF THE PHYSICAL PHENOMENON 

In order to quantitatively predict the physical phenomena one needs to describe 
it in mathematical terms and measurable physical quantities. Once a mathematical 
description is established (usually in the form of governing equations for the dependent 
variables) the solution of these equations are sought. An important intermediate step 
between the description and solution is the nondimensionalization of the governing 
equations. This not only simplifies the equation and in some instances even reduce the 
number of dependent variables, but also filters the natural physical parameters effecting 
the phenomena. The following sections concern themselves with the above issues with 
particular emphasis on the fundamental assumptions and approximations. 

3.1 Governing Equations 

Most real life fluid flow phenomena are mathematically represented by the well 
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known Navier-Stokes (N-S) equations which are based on the continuum hypothesis. The 
N-S equations are a set of nonlinear partial differential equations arrived at by the 
conservation of transport properties such as mass,momentum and energy for an 
infinitesimal control volume. In vector notation they are as follows: 

Conservation of mass, 

$2- + V.(pu) = 0 
01 - 


Conservation of momentum, 

p — + p(u*V)w = -Vp + V-t +/ 
dt 


( 3 . 2 ) 


Conservation of Energy in terms of the enthalpy. 


p— + p(u»V)ft = -V m q + 
dt 


dp 

dt 


+ (u«V)p 


;+ 4 > 


where 

. : is the divergence operator 
V : represents the gradient vector operator 
p : is the density 
p : is the thermodynamic pressure 


( 3 . 3 ) 


H : stands for the velocity vector 
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? : is the stress tensor 
7 : stand for additional body forces 

J 8 

h : is the specific enthalpy of the fluid 
q : is the heat flux vector 
$ : is the dissipation function 
The dissipation function is defined as, 

4> = ~ t • def( u ) ^ 

Here def(u) stands for the rate of deformation tensor and is defined as, 

defiu) = |[W+(V«r] (35) 

The superscript denotes the transpose of a tensorial quantity. Some points to 
note in the above equations are when the gradient operator ( v ) acts on a vector quantity 
it results in a tensorial quantity. The continuity and energy equations are scalar whereas 
the conservation of momentum is a vector equation out of which follow three scalar 
equations (assuming Euclidian space) depending upon the choice of coordinate system. 
The above set of equations together with the boundary conditions are necessary to solve 
the problem completely, but they are insufficient as there are more unknowns than the 
number of equations. The following section wrestles with this problem by invoking some 
fundamental assumptions and constitutive relations. 

3.2 Constitutive Relations And Fundamental Assumptions 
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The flow field is determined by the velocity ft and the thermal field by the 
temperature 7 , therefore constitutive relations are used to reduce the equations (3.1 )-(3.3) 
in terms of these dependent variables. Firstly the stress tensor ? is expressed in terms 
of the velocity gradients by assuming the working fluid to be Newtonian which along 
with the Stoke' s hypothesis is given by: 

t = \idef(u) - -^p/(V«w) (3.6) 

Here " / " denotes the unit tensor which is like the Kronecker delta in cartesian 

tensor notation. Substituting the definition of rate of deformation tensor def(u) * n 
stress tensor and then taking its divergence the conservation equation of momentum 
reduces to what is normally called the Navier-Stokes equations: 

p — — + p ( u*V ) w = -Vp + V-[^(W)]+V-[p(W)*]+/ -|v[»i(V*i?)] 0-1) 

dt 3 

Since the present study is not concerned with buoyancy effects the body forces due 

to gravity is assumed to be negligible and defining the pressure as: 

P = p + -^p(V»u) (3-8) 

equation ( 3 . 2 ) can be reformulated in the following form: 

p— + p ( m • V) u = -VF + V* [p ( W)] + V*[p (W)*] (3.9) 


As far as the conservation of energy equation (eq. 3.3) is concerned the first 
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simplification is done by using the Fourier heat conduction law which gives a relationship 
between the heat flux vector ( q ) and the temperature gradient is: 

q = - k [ V( T) ] ( 3 - 10 ) 

The specific enthalpy of the fluid h is changed to temperature by using the 
thermodynamic identity: 


Dh DT nr . Dp 

P = pc, + (1 - P T ) — ^ 

Dt p Dt Dt 


( 3 . 11 ) 


where 2 stands for the substantial derivative given by: 
Dt 


D ( ) = an 

Dt dt 


+ («-V)() 


(3.12) 


and p denotes the bulk expansion coefficient ( or thermal expansion coefficient) defined 

as: 


P ■ -- 

P 


a P 


8T 


(3.13) 


And it is zero for incompressible fluids and \]j for ideal gases. Substituting equations 
(3.11) and (3.13) in the conservation of energy equation (3.3) and with some algebraic 
manipulations it reduces to: 


pcff*pc p (u‘V)T = V.[tV(D]*pT 


dt 


+ $> (3.14) 


For an ideal gas the energy equation takes the form of: 
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P c p j^ + pc p (u-V) T = V.[*V(T)] + 


~~ + (w*V) J p 
Cf 




( 3 . 15 ) 


Equations (3.1), (3.7) and (3.15) are sufficient to solve for the flow and thermal 
field, but there are situations where they can be further simplified by the type of flow. 
In this study two such classes of flows have been studied, namely incompressible flow's 
and thermally expandable flows. The necessary condition for the establishment of 
incompressible flows is that the Mach Number is much less than one ( M < I) and if the 
temperature gradients are low in the domain considered the fluid properties ( p , p, P,A,c p ) 


are constant. Under this situation the continuity and momentum equations are greatly 


simplified with lot of terms dropping out such as dp/dt in e fl- (3.1) and v» [p( Viz )* ] 

in eq. (3.9). An important implication of this is that the continuity and momentum 
equations are decoupled from the energy equation, hence can be solved independently for 
the velocity field without worrying about the thermal effects. 

When there exists substantial temperature gradients in the domain at low Mach 
numbers then the fluid properties are no longer constant and vary with temperature (only). 
Such low speed compressible flows are called thermally expandable flows or anelastic 
flows. It usually implies the density of the working fluid varies only as a result of 
isobaric thermal expansion; in effect removing any acoustic phenomena from theoretical 
considerations. Under this situation no terms in equations (3.1) and (3.7) drop out unlike 
in the incompressible flow situation and further the energy equation now is coupled with 
momentum equations now due to temperature dependent fluid properties. 

If the Mach number is low ( M«1 ) an interesting formulation for the energy 
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equation (eq. 3.15) results, the second term on the right hand side of eq. (3.15) and the 
viscous dissipation function ($) becomes negligible and they drop out of the equation. 
Then the energy equation is reduced to: 

pcy|^ + p c/«-V)T = V-[kV(T)] < 316) 

The above equation states that the energy equation in low speed flows is strictly 
a balance between the convective and diffusive (conduction) processes. It should be 
noted that the above formulation is the same for both incompressible flow and anelastic 
or thermally expandable flow, but with the difference that in the former type of flow the 
energv equation is uncoupled from the other the conservation equations, whereas in the 
latter type of flow' all the equations are coupled together and therefore must be solved 
simultaneously. In the Stirling engine heat exchangers the flow speeds are very low- 
compared to the speed of sound (i.e low Mach number) and the above formulation for the 
energy equation suffices. 

Equations (3.1), (3.7) and (3.15) provide four equations for the four dependent 
variables U>V,P and T for a two dimensional cartesian or axisymmetric coordinate 
system. These equations are listed in the Appendix A in one generalized formulation. 
The complete set of partial differential equations as described in Appendix B together 
with the boundary conditions are necessary and sufficient to describe the fluid flow and 
heat transfer in oscillatory flows. The basic assumptions used to derive the theoretical 
equations are summarized as follows: 


- The fluid is a continuum. 
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Newtonian Fluid. 

Stoke’s hypothesis. 

Low Mach number. 

No body forces or gravitational effects. 

Axisymmetric or two dimensional geometry. 

Fourier heat conduction law. 

No internal heat sources. - No radiation heat transfer. 
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CHAPTER IV 

NONDIMENSIONAL PARAMETERS AND BOUNDARY CONDITIONS 

Nondimensionalization of the physical problem is one of the most important steps 
in the solution process of the governing equations. The nondimensional variables which 
arise from this process not only simplify the problem but also serve to reveal key physical 
aspects of the phenomena. Also normalization of the physical problem provides the 
natural scales for the problem as dictated by the boundary conditions, physical constants, 
and geometry. A dimensional analysis of the governing equations has not been presented 
in this report although the various nondimensional variables are described in detail along 
with their physical significance. 

As mentioned earlier (Chapter II) in this report there seems to be no consensus on 
the standardization of the dimensionless frequency in oscillatory flows and this is highly 
desirable in order to correctly interpret the results and for future reference. It is the aim 
of this chapter to clearly group the nondimensional parameters for oscillatory flows found 
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in the literature and physically interpret them so as to make a case for consistent usage 


in the future. 


4.1 v alensi Number ( y a ) 

There are the three different nomenclatures one comes across in the literature for 
the dimensionless frequency in oscillatory flows. This nondimensional variable is the 
natural outcome of oscillatory flows equations due to the presence of the unsteady term 
in the Navier-Stokes equations. To be precise, it weighs the strength of the time 
derivative term in the governing equations just as the Reynolds numbers weighs the 
relative strengths of the convection and diffusion terms. The three definitions of the 
above variable are as follows 

Valensi number, 


Va « 


D. 2 

» (y> 


Kinetic Reynolds number. 




D. 2 

® (y) 


Womersely number, 


T N 


Q) 


( 4 . 1 ) 


( 4 . 2 ) 


v 


( 4 . 3 ) 
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It should be mentioned that Kurzweg (1985a) with whose results the present 
numerical simulations are compared with has used DJA in the definition of Womersely 

number instead of D n in his investigations of oscillatory flows in a parallel plate 

h ' 

channel. 

From the above definitions of the Valensi number (Ka) an d Kinetic Reynolds 
number ( R e ) it is shown that they are the same but differ in the vocabulary. The 

motivation to define it as the Kinetic Reynolds number ( Re w ) can be traced to its 

similarity in structure with the well known definition of Reynolds number. The Valensi 
number can be physically interpreted as the ratio of viscous diffusion time scale 

( dJ/Av ) to the oscillation period ( i/ w ). When the Va is low implies l/w-co or the 

viscous diffusion is fast relative to the oscillation frequency and the velocity profile 
approaches the familiar parabolic shape as seen in steady flows. For higher y a due to 
high frequency the viscous effects do not have the time to diffuse across the duct before 
the convected fluid arrives, hence instead of the parabolic profile one sees the presence 
of a small Stokes layer near the wall where the viscous effects are concentrated. 

The Womersely number (a) on the other hand gives a very geometrical 

interpretation, it is the ratio of the width of the duct ( DJI ) t0 the viscous penetration 


depth y v /o) • 
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Despite the various meaningful definitions for the dimensional frequency, the 
definition of Valensi Number ( y a ) has been adopted by the Stirling Engine researchers 
and this definition will be used throughout this report. 

4.2 Maximum Reynolds Number (i?^) 

Maximum Reynolds Number is the second similarity parameter or the 
dimensionless variable which arises out the normalization of the governing equations. For 
oscillating flows it is defined as follows: 

r € = ( 4 . 4 ) 

“ p 

The velocity is scaled by the maximum velocity amplitude ( U ) instead of the 

mean velocity since the mean velocity is zero in oscillatory flows. 

The Maximum Reynolds Number ( r € ) and Valensi Number ( y a ) together make 

max 

up the dynamic similarity parameters for two dimensional oscillatory flows. 

4.3 Geometric Similarity Parameters 

This similarity parameters arise when the length scales are normalized as dictated 
by the geometry of the problem. Since this work concerns itself with plane flows or two 
dimensional geometries. The axial length (£,) of the Heat Exchangers are normalized 

with the Hydraulic diameter (£) ). Therefore the dimensionless axial length which reveals 

h 


itself after the normalization is h\D } ■ 
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4.4 Derived Similarity Or Nondimensional Parameters 

In oscillatory flows it is not unusual to find in the literature additional 
nondimensional variables, derived from the fundamental variables described above. 

4.4.1 Relative Amplitude of Fluid Displacement ) 

Relative amplitude of fluid displacement (y| ) is one such derived parameter and 

it is defined as the maximum axial fluid displacement during one cycle divided by the 
length of the duct. Essentially, it states how far the fluid is pushed into the duct 
compared to the axial length of the duct for one oscillation period (if the fluid oscillated 
inviscidly). 



time during the cycle. 
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From a different perspective A r indicates the volume of fluid displaced in one-half 

cycle divided by the total volume of fluid contained in the duct. Since it is an geometric 
similarity parameter it plays an important role in characterizing entrance effects. 

4.4.2 Strouhal Number ( Str ) 

Strouhal number is another derived parameter widely used in external flows as a 
nondimensional parameter for the frequency of vortex shedding. Analogously for internal 
oscillatory flows the frequency of fluid oscillation w is scaled by the U wn /D h to arrive 

at Strouhal number ( Str )• 


Str 


w 


(%> 

D. 


(4.6) 


By a little algebraic manipulations it can be shown that the Strouhal number is not an 
independent dimensionless parameter and relates the Valensi number ( y a ) and Maximum 
Reynolds number (R e ) by the equation, 

QU 


Str = 


4 Va 
/fe 


(4.7) 


Earlier Stirling Engine Heat exchangers were designed and operated based on the Str 


and R e values. 

max 
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4.5 Boundary Conditions For The Governing Equations 

The governing equations presented above mathematically describe a whole class 
of fluid flow problems, the way problems are differentiated from one another are by the 
application of boundary conditions. Since the equations are solved in a finite domain the 
boundary conditions for each of dependent variable needs to be specified along the 
boundaries of the domain. For the present situation four different types of boundary 
conditions were used to close the problem. All the boundary conditions for the particular 
boundary is based on the Figure 4. 1 . 

4.5.1 Solid Walls 

Along the walls by virtue of no-slip, the fluid velocity is zero assuming the walls 
are impermeable. For the energy equation the walls can either be maintained at constant 
temperature or be a source of constant heat flux. Mathematically this can be expressed 
as, 

^ = ^<r 0. T-- OR (4.8) 

where the partial derivative w.r.t n implies gradient normal to the wall. 

4.5.2 Symmetry Planes 

On symmetry planes or axis of symmetry the normal gradient of the tangential 
velocity and the normal velocity are zero as far as the momentum equations are 
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concerned. By symmetry the normal gradient of the temperature is also zero for the 


energy equation, physically this implies the symmetry plane behaves as an adiabatic wall. 


These can be formulated for the plane shown in figure 4.1 as: 

w = — = 0 (4.9) 

dy dy 


4.5.3 Inlet Plane 

On the inlet planes mathematically the boundary conditions are of Dirichlet type 
i.e. values of the dependent variables are specified. Thus 

U = U Sin(<x>t ) 

in mw v ; 

V =0 (4.10) 

in 

T=T in 

It should be noted that the above equations are based on the inlet plane as shown 
in Fig. 4.1. The U velocity is time dependent and varying sinusoidally whereas the 
temperature is fixed w.r.t time. 

4.5.4 Outlet Plane 

At the outlet plane one does not know the boundary conditions apriori the normal 
practice is to keep the domain long enough such that the diffusive fluxes are negligible 
normal to the plane. Additional physical constraint for the momentum equations is 
derived by observing that for incompressible flows the mass fluxes are conserved. The 
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mathematical formulation for the boundary conditions are: 

W = = 0 (4.1 1) 

dx dx dx 

For oscillating flows the inlet and outlet plane are reversed after each half cycle 
due to the zero mean flow restriction. For instance if the cycle begins with the flow 
entering from the East side after half period the flow enters from the west side. 

4.5.5 Solid-Fluid Interface 

This boundary' condition is needed in the conjugate heat transfer problem, where 
there is an interaction between the fluid within the channel and the surrounding solid 
region. At a solid-fluid interface the heat flux across the interface is conserved by energy 
conservation principle and the temperature is continuous. These two boundary conditions 
have been implicitly implemented in the code since the solid-fluid domain has been 


solved together. 
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CHAPTER V 

NUMERICAL SOLUTION TECHNIQUE 

The governing partial differential equations (PDEs) which are generally 
elliptic in nature, are not tractable to analytical solving procedures. They are numerically 
integrated by one of the many discretization procedures. In this chapter,the solution 
methodology is discussed. 

Currently, there are numerous methods and ways to solve the partial 
differential equations arising in fluid mechanics, some of the existing popular and fairly 
standard methods are based on one of the following discretization schemes: 

i) Finite Volume Methods (FVM) 

ii) Finite Element Methods (FEM) 

iii) Finite Difference Methods (FDM) 

The underlying principle behind all these methods are one and the same, 
which is dividing the domain into smaller subdomains and thereby reduce the partial 
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differential equations into a set of linearized algebraic equations. The differences among 
the methods arise from the way these algebraic equations are arrived at from the original 
PDEs. Since no single method has till now been proven to be superior to the other, the 
finite volume method (FVM) has been chosen in this study. One of the distinct 
advantages of the finite volume method is that it lends itself naturally to the PDEs arising 
in fluid mechanics problems. 

The code developed to solve the governing equations is based on the 
research code called C.A.S.T (Computer Aided Simulation of Turbulent Flows) developed 
by Peric and Scheuerer ( See Peric et al. ). The original code limited in its ability to 
handle variety of thermal and time dependent boundary conditions has been broadened 
to include these type of flows. Further the numerical formulation of the energy equation 
has been revised to handle conjugate heat transfer problem such as that occurring in the 
regenerator of the Stirling engine. The objective of the following description is (i) to 
briefly describe the discretization technique (ii) solution technique employed to solve the 
algebraic equations and (iii) discuss the convergence criterion employed. 

5.1 Principle Of Finite Volume Method 

The finite volume method in general is based on the conservative property 
of the partial differential equations (PDEs) since the equations which themselves are 
derived from the conservation of certain physical quantities. This important attribute of 
the equations makes it possible to collapse all the individual equations into a generalized 
transport equation, thus facilitating one common algorithm for all the PDEs. Therefore 
for any generalized scalar variable <j)the transport equation can be written as: 
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d(P<t>) 

dt 


a*, 


i a 


a*, _ 




( P (/<f>-iy^) + — ^(pr-^-r^) = 

ox r " or v dr 


(5.1) 


Here the value of n determines the type of axes system used, when 0 
then the transport equation reduces to the cartesian coordinate system with the 
independent variable r being changed to the more familiar y. And, when 1 the 
equation represents the transport equation for a axisymmetric coordinate system. All the 
individual differential equations for each conservation equation (e.g. mass, momentum) can 
be recovered from equation (5.1 ) by choosing appropriate physical quantity for <k r and s 


which are given in Table 5.1 below: 

Table 5.1: Interpretation of (h ,r and 5 in the transport eq. (5.1) 

4> <p 


Equation 

4> 

r * 

h 

Continuity 

1 

0 

0 

x-Momentum 

u 

M 

3P 3, dU. 13. „3K 

r-Momentum 

V 

M 

dr dx dr r n dr ^ dr r 2 ' 

Energy 


k 

0 
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The general discretizing procedure for the transport equation given in eq. 
(5.1) is as follows. In the first step, the computational domain is divided into small 
rectangular control volumes. The grid points where the dependent variables are solved 
for lie at the geometric centre of the control volume. Then the transport equation 5.1 is 
integrated over each control volume. By applying Gauss’s theorem to the integral 
equation results in a integro-differential equation. Physically the integro-differential 
equation is a relation between the net increase of the considered quantity per unit time, 
the total net convective and diffusive fluxes across the control volume boundaries and the 
source (or sink) terms within the control volume. This integro-differential equation is 
then discretized with some assumptions and linearization to arrive at the linear algebraic 
equation for the control volume. If the above process is repeated for all the control 
volumes one arrives at a set of linear algebraic equations for the transport equation 5.1 . 
The natural appearance of the fluxes at the control volume faces makes the whole scheme 
globally conservative. These set of equations are then solved to give the values of the 
considered quantity at the grid point locations. Special coupling algorithm is used to link 
all the sets of algebraic equations for a given physical quantity. A brief outline of the 
above described is as follows for a detailed explanation the reader is referred to Peri6 and 
Scheuerer (1989). 

The outline will be presented with reference to the Cartesian coordinate 
system, which means in the transport equation (5.1) n“0 and r=v . 



INI.NJI 


3a 




Figure 5.1. Control Volume arrangement and grid numbering. 
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5.2 Control Volume Variable Storage 

Figure 5.1 show a part of the integration domain which is subdivided into 
small rectangular control volumes by the intersection of x and v lines. The code used 
in this study uses the so called co-located variable grid arrangement. In the co-Iocated 
variable arrangement all the dependent variables ( u , y , 7 ;.. .etc) are stored in the same 

location as shown in figure 5.2. Patankar (1980) had shown how the co-located variable 
arrangement gives rise to oscillatory (checkerboard) solutions but a special interpolation 
procedure is used in the present code to determine the cell-face velocities (the main 
reason for unphysical solutions) to suppress the checkerboard solution. In the discretiza- 
tion scheme to be discussed next, a typical control volume (CV) containing the grid point 
P (see figure 5.2) is integrated. All the surrounding grid points are identified by their 
sense of direction relative to the grid point P, such as the grid point E located to the right 
of point P. All the quantities calculated at the CV faces are denoted by lower case 
subscripts such as "e" for the quantity calculated at the east side cell face. The open 
arrows denote the mass fluxes at the CV faces in the x and v directions. 

5.3 Integro-Differential Equation 

The first step in the discretization of the transport equation (eq. 5.1 ) is to 
integrate it over a control volume ( 5 V ) to yield: 


f dv+f [— (p (/ 4 ,-r ' .§£) 

Jfiv Q( J&v Qx * dx 




dy 


(pVi> 


r.— )]<Jv» f 

* dy Jt 


6 V 


[S,]dv (5.2) 


Where 
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6v = £i.x{jy { 53 ) 

is the two dimensional "volume" element in the cartesian coordinate system. The second 
term in eq. (5.2) is the volume integral over the divergence of the convective and 
diffusive fluxes. Applying Gauss’s divergence theorem it is transformed into a surface 
integral resulting in: 

f 6{ ^ dv + f "[ (p -T.— ) - (p £/<(> -T.— ) ]dy 

Jcv & Js * *dx *dx' v 

* O (p ^ - r »f - ( p y * - r *f : U* ,5 ' 41 


■ Ow 

A closer look at equation (5.4) reveals that it is nothing but a balance between the 
rate of accumulation of <j) within the CV and the net transport of (J> by convection and 

diffusion across the CV faces plus the source or sink terms within the CV. Further the 
eq. (5.4) is still exact in the sense no approximations have yet been introduced. The next 
step is the discretization of the eq. (5.4) thus introducing the approximations. 

5.4 Discretization Scheme 

The discretization of eq. (5.4) is done in two steps. In the first step the surface 
and volume integrals appearing in are approximated by utilizing the mean value theorem. 
The linearization of the coefficient are also done in this step. In the second step the mean 
values of various transport quantities arrived at in the first step are discretized and related 
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to the CV grid points. 

5.4.1 First Level Approximations 

The approximations for the integrals appearing in eq. (5.4) using the mean value 
theorem involves two crucial assumptions, viz: 

► the fluxes through the CV faces are approximated as the product of the mean flux 
per unit area going through the centre of the particular C V face and the CV face 
area. 

► the integrals involving the time derivative and the source terms are approximated 
as the product of the mean value of the integrand associated with the CV centre 
and the CV. 

► the coefficients are also linearized in this step. All the mass fluxes and the 
diffusivities are evaluated with values from previous iteration. 

Using the above three assumptions eq. (5.4) can be written as: 


d(p‘<t>) 

dt 


tv * - [/„ 




dx 


<t> -T\ (-^) by] 






5 


(5.5) 


= V v 

Where 


m r = pU(by) ; rn y = pF(8x) (5.6) 

are the mass fluxes across the control volume faces in the x and y directions 
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respectively. Also the starred ("*") variables indicate values obtained from previous 
iteration. 

The second step of approximations involves the discretization of the convective 

fluxes, diffusive fluxes, source terms and time derivative. 

5.4.2 Discretization of Convective Fluxes 

This step involves the discretization of the convective fluxes or the mass fluxes 

found in eq. (5.5). Special care has to be taken to find the values of the dependent 

variables at the CV face (e.g. A , etc) with respect to the values at the grid nodes ( i.e 

* € 

the centre of the CV). The code used in this study uses a flux-blending approach or a 
hybrid method to evaluate the variable at CV faces which is linear combination of two 
methods with different order of accuracy, namely 

► the first-order Upwind Differencing Scheme (UDS), and 

► the second-order Central Differencing Scheme (CDS). 

A pure CDS cannot be used in flow situations where there is a flow reversal and 
recirculation as it can give rise to unphysical solutions (see Patankar, 1980), whereas UDS 
albeit lower in accuracy ensures a diagonally dominant, positive definite coefficient 
matrix. This is achieved in UDS by replacing the cell face value by the grid point or 
node value closest to the CV face depending on the flow direction, unlike in CDS where 
the cell face value is evaluated from a linear interpolation of adjacent grid point values. 
Using the hybrid approach the value of the dependent variable at the east cell face can 


be written as: 
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♦. - i> e VDS + Y( < t> € C£,S-< t ) e l/D5 )* < 5J > 

Here y weighs the contributions of UDS and CDS and is 0 for a pure UDS and 1 for 
pure CDS. The implies the values evaluated at a previous iteration level thus making 
it a "deferred correction" approach. The deferred correction ensures the diagonal 
dominance required for the solution of the algebraic equations even for a pure CDS. In 
the present study a pure UDS ( y = 0) was used for all the cases investigated. 

5.4.3 Discretization of the Diffusive fluxes 

The diffusive fluxes are the first derivative terms (gradient terms) in eq. (5.5) 
multiplied by the diffusivity (r ). A second order central difference discretization 

method is used to evaluate these terms and a typical form of it for the gradient at the 
"east" cell face would be: 


d<K _ 

6x e Ax. 


(5.8) 


It is important to note the <£ are evaluated at the grid node value (upper case). Similar 
formulations can be derived for the gradients at other cell faces. 

5.4.4 Discretization of the Source Terms 

The source terms appearing in the right hand side of eq. (5.5) is replaced by the 
value obtained at the centre of control volume (or grid node P, see fig. 5.2), i.e.: 


*V v c V 6v 


(5.9) 


In case of non-linear special care has to be taken to linearize the source terms such that 
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only the positive contributions are added to the algebraic equation. For further details 
regarding this see Peric and Scheuerer(1989). 

5.4.5 Discretization of the Time Derivative 

Two basic assumptions are utilized to discretize the time derivative term in eq. 

(5.5). First, CV mean value of 4 , is approximated as the grid node value and 

secondly, the value of dependent variable is assumed to be varying linearly with time. 
Utilizing these two assumptions one gets: 

6 (p<t>) P4> r - P°p<> )0 p (5.10) 

dt “ t - t° 

The superscripts "0" stands for the values obtained from previous time step. 

Since all the space derivatives and source terms are evaluated at the new or 
current time level it makes the whole discretizing scheme fully implicit. In other words 
there is no restriction on the time-step chosen. 

5.4.6 The Final Form of the Discretization Equation 

After substituting for the approximations described in Sections 5.4.1-5.4.6 into the 
eq. (5.5) one gets the Final discretization equation for grid node P as: 

<Vi> P = c«4v + a &E + a s4>s + a r&N + b p (5II) 

A look at eq. (5.11) reveals that the value of the dependent variable at the grid node P 
is dependent on the surrounding grid nodes. The coefficients "a" contain the 
contributions from the convective and diffusive fluxes and " contains the source term. 
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The detailed formulations of the coefficients can be found in Perfc and Scheuerer (1989). 

5.5 Solution Algorithm 

The governing equations presented in Chapter m (eqs. 3. 1-3.4) are all coupled 
together for a general fluid flow problem. For incompressible low speed flows the 
momentum and continuity equations are strongly coupled together. Hence special 
algorithms need to be used to solve for the dependent variables, in the numerical code 
used in the present study the well known and tested SIMPLE algorithm (Semi-Implicit 
Method for Pressure-Linked Equations) has been implemented. A brief outline of the 
SIMPLE algorithm will be presented next. For unsteady flows the solution algorithm is 
applied to each time step. 

5.5.1 SIMPLE algorithm 

The SIMPLE algorithm is sequential step by step solution procedure where each 
of the governing equations (e.g. continuity ,x-momentum, etc.) are solved one after another 
and then coupled together by physically derived algebraic relations. The algorithm 
consists of the following steps: 

0) First step consists of initialization of all the dependent variables such that the 
finite volume coefficients (the fluxes) and the pressure difference (source term) 
in the momentum equations can be evaluated. Any sensible initial guess value can 
be used, for unsteady flows such as the present problem values from previous 
time step can be used as a good initial estimate. 

1) Next the finite volume coefficients of the x-momentum equations are assembled. 
Then the resulting set of linear algebraic equations are solved to yield the axial 
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velocity jj* . Since the resulting new velocity is only approximate (based on 

initial guessed pressure), the algebraic equations have been solved by an iterative 
solution algorithm instead of direct matrix inversion algorithm. These iterations 
are called "inner" iterations and the number of these "inner" iterations can be user 
controlled. 

2) The same procedure (Step 1) is used to obtain the normal velocity y from the 
y-momentum equation. 

3) Since the initial pressure is guessed, the velocities (/* and y* will not satisfy 

the continuity equation even though they will satisfy the momentum equations. 
In this step a pressure correction equation is derived (see Appendix B) to estimate 
the pressure and its associated velocity and mass fluxes from the continuity and 
momentum equations. These corrected velocities will then satisfy the continuity 
equation but will throw the momentum equations out of balance in the process. 

4) Once the velocities satisfying the continuity equations are found all the scalar 
transport equations are solved in the present case the energy equation. The 
energy equation is solved in the same process as described in Step 1. As the 
energy equation is decoupled only one "inner" iteration is performed, although t 

could have been solved after the true velocities have been found. 

5) The residual norms are computed for all the conservation equations which in an 
ideal case of correct solution should go to zero. These residual norms are 
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normalized by appropriate reference quantities. If any of the normalized residual 
norms is greater than the user specified convergence criterion the algorithm 
returns to Step 1 and uses the current values to evaluate the new finite volume 
coefficients. 

6) If all the normalized residuals norms are smaller than the specified convergence 
criterion than convergence is declared. For unsteady problems, the time counter 
is incremented by the time step size ( ) and the algorithm return to Step 0 with 

the initial guesses for the dependent variables taken from the previous time step. 
In the present study the convergence criterion was choosen to be 0.1% of the 
reference residual norms. Further, special care has to be taken for oscillatory flows when 
evaluating the reference mass flux for the residual normalization(Step 6). Since at the 
instant of flow reversal the mean flow velocity is zero. In these situations the code has 
been modified to use the reference mass flux based on the maximum inlet velocity 

( U )• 

max 

5.6 Code Modifications For Oscillating Flows 

The computer code CAST has been modified to account for the cyclical nature of 
the oscillatory flows and the switch in the boundary conditions. Although care has been 
taken to retain the structure of the code and utilize the Vectorizing capabilities of the 
Cray YMP supercomputer. Briefly, few of the major changes include: 

► At flow reversal (0°, 180°, 360°) the inlet and exit boundary planes are switched b 
account for the zero mean flow situation. 

► A new energy equation (assembly & evaluation of the FV coefficients) routine has 
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been written along with the proper boundary condition to solve for the heat 
transfer problem. 

- The code has also been modified to solve the conjugate heat transfer. The change 
has been made such that the numerical code is transparent to the presence of solid 
and fluid region (i.e. the energy equations is solved together for the solid and 
fluid regions. This has been achieved through modifications of the diffusive 
coefficients and source terms details of which can be found in the book by 
Patankar (1980). 

In addition to these, minor modifications have been done to accelerate the 

convergence of the equations. 
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CHAPTER VI 

OPERATING CONDITIONS AND GEOMETRIC MODELING 

This chapter concerns itself with the nondimensional parameters conditions under 
which the Stirling engine heat exchangers operates. This chapter also addresses the way 
the heat exchangers are modelled numerically such that the geometry closely 
approximates to the one found in the engine. The operating parameters are given 
specifically for the NASA SPRE (Space Power Research Engine), the object of present 
study. Figure 1.2 showed how the heat exchangers in the SPRE are located, the cooler 
starts from the compression space and opens into the regenerator which in turn is 
connected to the heater and which opens into the expansion space. 

6.1 Operating conditions 

The parameters under which the SPRE operates were taken from the one 
dimensional code GUMPS, which simulates a Stirling cycle engine and has been used 
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extensively by researchers at NASA. The numerical simulation of the heat exchangers 
used in the present study were based on GLIMPS and these are listed below in columar 
format for the Regenerator, Cooler and Heater respectively. 

Regenerator 


PARAMETER 

Hyd. Diameter D h , (m) 

GLIMPS CODE 

0.133 3x10' 

PRESENT STUDY 

Wire Diameter D w , (m) 

2. 54x1 O' 5 

5.0x1 O' 3 

Matrix Length L, ^ (m) 

24.63x1 O' 3 

12.00 

Mean Pressure (Mpa) 

f J 

15 

n/a 

Gas Temp. (Hot side) (TO 

617.2+(12.6) Sin (to t+2 . 

.07) 274 

Gas Temp. (Cold side) (^C) 

345.1 +(2.10) Sin (to t+2 

.31) 273 

Valensi Number (Va) 

1.5 

^5^00 


265 

^75,12000^ 


0.478 

0.25 

Cooler 

PARAMETER 

GLIMPS CODE 

PRESENT STUDY 

Hyd. Diameter D h , (m) 

1.524x10 3 

5.0x1 0^ 

Tube Length L, (m) 

95. 2 5x1 O' 3 

6.0 

Number of Tubes & 

1584 

1 

1 

Mean Pressure (M^a) 

15 

n/a 
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345.1 +(2.1) Sin (tot-2 .31) 340,350 

332.7+(12.9)Siu(wt+0.13) 330,340 

324.3+7.5xl0 5 Si.n (<»> t+0 .7" 

350 

300000 15000 ,30000- 

& 60000 

0.686 0.714 



Heater 



PARAMETER 

GLIMPS CODE PRESENT STUDY 

Hyd. Diameter D h , (m) 

1.27x10 3 

5.0x10 3 

Tube Length L, (m) 

90. 17x1 0' 3 

7.0 

Number of Tubes ^ 

1632 

1 

Mean Pressure (M^a) , 

15 

n/a 

Gas Temp. (Hot side) (TK) 

625.8+(29.1) Sin (to t+0 . 33 ) 

630,620 

Gas Temp. (Cold side) (4c) 

617.2+(12.6) Sinito t+2 . 07 ) 

610,620 

Wall Temp. (°K) 

324.3+ 1.4x1 O' 4 Sin (to t+0 . 76 ) 

650 

Valensi Number ( Va ) 

88 

44,88 & 176 


16500 

8250, 16500 



& 33000 

^x 

0.686 

0.714 
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The regenerator parameters used in this study are diffrent from that predicted by 
the GLIMPS code since the regenrator was modelled as a seperate entity. Whereas the 
operating conditions for the regenerator got from GLIMPS is based on a complete 
simulation of the Stirling engine i.e. including the cooler and heater. Also the wall 
temperature fluctuations were assumed to be zero in the present study as they are small. 

Working Fluid 

The working fluid used in the numerical simulations was helium which is the same 
as used in the SPRE. The following properties of helium at standar atmospheric pressure 
were used for all the numerical simulations conducted in this study: 

Density (kg/m 3 ) 0.200 

Dynamic Viscosity (N.s/m 2 ), 2.83x10 5 

Specific Heat, C p (J.kg’/TC 1 ) 5200.00 

Thermal Conductivity, ic,fW.m' l .' , K) 0.20439 

Prandtl Number ( Pr ) 0.72 

Regenerator Metal 

The regenerator metal used in the conjugate heat transfer problem was chosen to 
be aluminum and the following properties at 20°C were used in this study: 


Density (kg/m 3 ) 

/ 

2707.0 

Specific Heat, C p (J.kg / 

K- 1 ) j 

896.00 

Thermal Conductivity, jc,(W.m'7TC) 

20.4 


OHK3NM. 

nr POOR QUALiTY 


OF POOR 


54 


Table 6.1 : Test cases investigated in the present study for fluid flow analysis. 

























































Figure 6.1. Envelope in which different Stirling Engines operate, together with: i) 
Criterion for transition from laminar to turbulent flow, ii) Different test 
cases studied in the present work. 
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Table 6.1 lists all the cases studied in the present work. All the cooler cases begin 
with the letter "C", the heater with "H" and the regenerator with "R" respectively. Figure 
6.2 shows these cases on a R e m , v and Va plot with various transition criterion found in 
the literature for laminar to turbulent transition. As it can be seen the SPRE heater and 
cooler heat exchangers operate on the transition regime. But a laminar analyis is still 
pertinent since the flow inside the heat exchangers is part laminar and part turbulent 
within an oscillation cycle. 
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6.2 Geometric Modelin 


Figure 1.2 showed an conceptual geometry of the SPRE, the heater and cooler 
components consists of bundles or*mbes and the regenerator is a foil type matrix. The 
heater and cooler are modeled as circm^r tube with a finite length. Only one tube for 
each component was considered since a thrcmgjh analysis involving all the tubes is beyond 
the scope of the present study. The regeneratoRis one of the most important in the 
Stirling engine and is also one of the most difficult to^riodel since it is a matrix. In the 
present study this is resolved by modeling it as a parallel-^ates channel with geometric 

finit^t 


similarity i.e. the matrix is replaced by a solid plates of 
hydraulic diameter. In summary : 


thickness with the same 


Heat Exchanger 


n 

; i 

Geometry Mod^elhjd 


Regenerator 

Heater 

Cooler 


Parallel-Plates Channel 
Circuler Tube 
Circular Tube 
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CHAPTER VII 

OSCILLATING FLUID FLOW ANALYSIS 

In this chapter results for the effect of flow oscillations on the velocity field and 
associated friction losses are presented and discussed. The investigation has been carried 
out for a wide range of nondimensional parameters with an eye on the Stirling engine 
heat exchangers operating conditions. The cases investigated and their operating 
parameters are given in Table 6.1 (Chapter VI). 

Any numerical discretization method gives rise to so called truncation and 
discretization errors. Hence, code validation is an important element in any numerical 
simulation, in the present study the results of the numerical simulations have been 
compared with existing analytical and experimental efforts. The first two sections 
(7. 1,7.2) deal with this aspect concurrently with the fluid flow analysis. 
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7.1 Fluid Flow in Parallel-Plates Channel and Code Validation 

Since the regenerator is modelled as a parallel-plates channel the fluid flow results 
for this component are presented next. 

7.1.1 Analytical Solution For The Velocity 

Kurzweg (1985a) analytically solved the N-S equations governing the oscillatory 
for flow between two-parallel-plates channel. The geometry investigated is shown in 
Figure 7.1, where an array of parallel-plates channel are connected at ends, to reservoirs 
maintained at different temperatures (only part of the geometry as simulated numerically 
is shown). The flow is set to motion inside the channel by a sinusoidally varying 
pressure gradient and a temperature difference between the end reservoirs ensures that a 
constant axial linear temperature gradient is maintained along the channel throughout the 
cycle. Also, he assumed the walls of the channel to be thick thus signifying a conjugate 
heat transfer problem or a special case of the generalized constant heat flux thin wall 
problem. 

The channels were assumed to be long such that the fluid flow is "fully 
developed" or to be more precise, the axial velocity profile is constant along the channel. 
Under this assumption the momentum equations simplify considerably — in fact, only the 
axial momentum remains — and are tractable to analytical solution techniques. By 
neglecting the initial conditions he solved for the quasi-steady axial velocity distributions, 
i.e, the velocity distributions does not change from one cycle to another cycle at any 
instant in the cycle and is given by: 

Where y = is the imaginary unit, # denotes the real part of a complex 
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U{T],t) =U 0 f(r\) e iuf = u 0 st 


cosh ( y/Igr| 
cosh ( Jla ) 


quantity, jj q is an arbitrary velocity scale, rj = y / a the normalized coordinate distance 
[ y] normal to the flow direction, t the time, a = g/u/v the Womersley number 

or the nondimensional frequency (note va=4a 2 )> and X = Idp/dx^ a 2 / p^ 0 v ‘he 
nondimensional amplitude of the imposed sinusoidal pressure gradient. Once the velocity 
distribution is established practically useful quantities such as the wall shear stress can 
now be established. The temperature profile is discussed in the heat transfer section to 


be presented later in this report. 

7.1.2 Reformulation of The Analytical Solution 

In order to compare the velocity profiles by the numerical solutions w ith 
the analytical solution presented above special attention has to be given to carefully match 
the boundary conditions. A direct comparison with the above solution could not be made 
because of the different boundary conditions and the two-dimensionality of the numerical 
simulations. In the derivation of the analytical solution the parallel -plates channel was 
assumed to be infinetly long (thus rendering it a one dimensional) and more importantly 
the flow was established by applying a sinusoidally varying pressure gradient 
[dp /dx = | dp / dx\ MX cos ( w t )]• But in the numerical solution the domain is two 
dimensional and finite in length with the flow being established by a sinusoidally varying 
velocity [ p = p sin ( u) t) ] at the inlet, the reason for choosing a velocity boundary 
condition instead of a sinusoidally varying pressure boundary condition was the ease of 
implementation numerically. Hence the analytical solution for a periodic velocity 


6\ 



Figure 7 . 1 . The geometry for parallel-plates problem (conjugate heat transfer) with the 
inlet and outlet boundary conditions. 
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boundary condition has to be derived first and the proper axial location chosen. The axial 
center of the channel was chosen for comparing with the theory since it is far away from 
both the ends the effects of finite channel length can be assumed to be negligible at this 
location. The solution for the velocity profile [Eq. (7.1)], essentially, stays the same but 
a phase difference has to be added to account for the velocity bounadary condition. The 
following equations outline the procedure to derive the phase difference 4 ^ to be added 
in order to correspond to the numerical simulations: 

Given, 


U ln « u MX sin (<*t) 


( 7 . 2 ) 


Now defining the y given in Eq. (7.1) to be the same as cr max (maximum inlet velocity), 
the reference parameter used in numerical solution, or: 


Uo s tWx 


( 7 . 3 ) 


and with little algebraic manipulations the nondimensional pressure gradient \ in Eq. (7.1) 
can now be related to the known numerical parameter a (or ya ) by the expression: 


Here 5 is a complex quantity and the double vertical lines ( || || ) stands for the 
absolute value or modulus of the complex quantity. Its given by the expression, 




tanhi JTa 


(7 . 5 ) 


Also the phase difference 4 ^ to be added to the analytical solution for a sinusoidally 
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varying velocity input is given by: 


<t) u = 90 + 


tan 1 


*<»«> 

* (#«) 


( 7 . 6 ) 


Where symbols and denote the real and imaginary parts of a complex number. 
Finally, the velocity distribution corresponding the boundary condition given by Eq. (7.2) 


is given by: 


m n . t) « cu f a ^ - iu* * 


cosh(v'icr)) 
cosh (fia) 


i { w c * (f'j. * 


( 7 . 7 ) 


Where f (^) is a complex function varying with the normal distance ( v ) and given as: 


f (n) ■ ri - ( 7 . 8 ) 

11 ^ ||6 U « cosh(JIZ) 

Equation (7.7) was used to compare with the numerical velocity predictions at the mid- 
plane of the channel as the entrance effects due to the finiteness of the computational 
domain can be assumed to be negligible at this axial location. The practically useful 
quatities such as the wall shear stress ( Xfcr ) (directly related to the friction factor) and the 
pressure drop across the channel were also compared with the numerical predictions. The 
corresponding equations for these two quantities can be easily derived from their 
definitions and for a sinusoidal velocity input are given by: 


12 


( 7 . 9 ) 


j^<>xU m&x sin(wt+^) 

Here <j is called as the wall shear stress augmentation factor which can be shown to be: 
and <{> the lead phase angle to be added to the wall shear stress can be expressed as : 
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°X = 


z T = -i [ iji a - tanh (JTa) } 


\7 . iOa , 7 . iob; 


<t> T = <t> a " tan ~ X 


' % ( z t ) ' 

«M Z, ) 


(7.11) 


The instantaneous pressure gradient across the channel can be obtained by integrating the 
axial momentum equation across the channel and with little algebraic manipulations one 
obtains: 


w t + 4> p ) (7.12) 

\ V h ) 

Here 0 and <j) p are the pressure gradient augmentation factor and lead phase angle 
respectively, and are given as: 







[O t COS<t> T ] 2 + 


+o T sin< t>J 2 


(7.i3: 


4>p = tan' 1 


va 

12 


+ o t sin <|) T 


o T COS<J) T 


(7.14) 


Few interesting points need to be mentioned about the effects of Valensi number 
on the auge mentation factors and phase angles for the wall shear stress and pressure drop. 
When the Valensi number ( y a ) is low the augentation factors ( 0t and Cp ) approaches 
unity and reduce to the familiar steady state formulations. For high va the phase angles <j> 
and <f, p asymptotically reach 45° and 90° respectively or the wall shear stress and pressure 
drop are out of phase with the inlet velocity by these angles. In addition for high y a the 
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augmentaion factors <j^ and o p reach as high as 8 and 120 times the steady state value. 
Figures (7.2) and (7.3) graphically demonstrates these effects, the augmentation factors 
are show in Fig. 7.2 and the phase angles in Fig. 7.3 both of them plotted versus the 
Valensi number ( va )• 

Equations (7. 7), (7. 10) and (7.11) were used to compare the axial velocity profile, wall 
shear stress at the axial mid-plane and the instantaneous pressure drop across the channel 
for the cases R, and R 2 . The results of these comparisons are discussed next. 

7.1.3 Comparison with Numerical Simulation 

Axial Velocity 

The regenerator cases (R, and R 2 ) were chosen to validate the code against the 
analytical solution. Fig. (7.4) shows a plot of the normalized axial velocity ( [//[; ) 

vs. the nondimensional transverse distance (y/ a ) f° r the case Ri Re^* ^ anc * Va ”2.5 
at different velocity phase angles (from 30° to 360" with 30° increment). The symbols 
are used for the analytical solution (Eq. 7.7) and the dotted lines are used for the 
numerical predictions. It should be appropriate to mention that the comparison was made 
at the axial mid plane { x - L/ 2 ) as shown in Fig. 7.1. The velocity profile exhibits the 
familiar parabolic profile due to the low frequecy (Valensi number) and is completly in 
phase with the input velocity. The agreement between the analytical solution and the 
numerical prediction is excellent. Fig. 7.5 shows a similar plot for the case R 2 ,j?e “ 
12000 and y a = 400. Here one can observe the effect of the high \r a on the velocity 
profile, by the presence of a small Stoke’s layer near the wall and with the flow field 
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almost uniform in the channel core. The rationale behind this phenomena being that the 
flow reversal (switch in the flow direction) takes place before the viscous effects have 
time to diffuse down along the radius. Also it can be observed from the Figure (7.5) the 
flow reversal is also captured very accurately by the numerical code even for high y a . 
Again it can be seen that the agreement with the analytical solution is excellent. 

Wall shear stress and Pressure Drop 

Figures (7.6) and (7.7) shows the normalized wall shear stress versus the velocity 
phase angle at the axial mid plane for the cases R, and R : respectively. The 

normalization factor chosen was the coefficient of the sine function in Eq. (7.9). The 
symbols denote the analytical result (Eq. 7.9) the solid line represents the numerical 
predictions. When the Valensi number is low, the -j (wall shear stress) is in phase with 
the inlet velocity phase angle and its magnitude is exactly equal to the steady state value. 
But for high y a (case R 2 , Fig. 7.7) the t w is not only 45° out of phase with the inlet 
velocity but also its magnitude is augmented four times the steady state value. 

Figures 7.8 and 7.9 show similar figures for the normalized pressure drop along 
the channel plotted versus the velocity phase angle , again the symbols are for the 
analytical solution (Eq. 7.12) and the solid line for the numerical prediction. Fig. 7.8 is 
for the low valensi number case (R,) and Fig. 7.9 for the high y a case (Rj). But unlike 
the wall shear stress the pressure drop for high y a is 90° out of phase with the inlet 
velocity and its magnitude is 40 times the steady state value. 
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Figure 7.4. Comparison between analytical and numerical velocity 
profiles at different velocity phase angles for the conjugate heat transfer 
problem, Case R, .Re^-15 and Va-2.5. 

[Symbols : Analytical-Kurzweg ;; Dotted Lines : Numerical ] 
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Figure 7.5. Comparison between analytical and numerical velocity 
profiles at different velocity phase angles for the conjugate heat transfer 
problem, Case R : 2000 and Va-400. 

[Symbols : Analytical -Kurzweg ;; Dotted Lines : Numerical ] 
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Figure 7.6. Comparison between Analytical and Numerical normalized Wall 
shear stress ( x u ) versus the velocity phase angle for Case R. : Re m 75 
and Va ~ 2.5. 
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Figure 7.7. Comparison between Analytical and Numerical normalized Wall 
shear stress ( x w ) versus the velocity phase angle for Case R 2 : “ 12000 

and Va m 400. 
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Figure 7.8. Comparison between Analytical and Numerical normalized 
Pressure Drop (A P/ L) versus the velocity phase angle for Case R, • Re 
75 and Va - 2.5. 1 
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Figure 7.9. Comparison between Analytical and Numerical 
Pressure Drop (A P/L) versus the velocity phase angle for 
12000 and Va - 400. 
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7.2 Comparison with Experimental Data 

This section is concerned with the comparison of the numerically predicted 
solutions with the experimental work carried out at the University of Minnesota. For the 
sake of clarity a brief description of the oscillatory flow rig at the University of 
Minnesota is also included. 

7.2.1 Description Of the Experimental Setup 

The experimental study by the researchers at University of Minnesota was 
initiated to help understand the thermo- mechanical energy losses in the free piston Stirling 
engines in order to come out with better designs. Their preliminary survey suggested, a 
better understanding and characterization of the laminar to turbulent flow transition in 
oscillatory flows (see Sueme and Simon,1988). The initial efforts were concerned with 
understanding the mechanisms by which transition takes place and generally character- 
izing the fluid mechanics of oscillatory flows. Recently (see Seume et al.,1992), detailed 
measurements were carried out at a particular operating point, namely that of the heater 
tubes of NASA’s Space Power Research Engine(SPRE). The velocity measurements were 
taken at four axial stations located along the test section as shown in Fig. (7.10). Figure 
(7.10) also demonstrates the flow oscillation in the test section 9or tube) being affected 
by the piston-rod assembly. The test section is a circular pipe connected by smooth 
nozzles at both ends; the smooth contour of the nozzles ensures no flow separation upon 
entry. One end of the test section is connected to the flow delivery section the other end 
opens out to the room. The inflow conditions at the two ends are nearly symmetric and 
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the axial locations are represented by s / cf where s is the distance from the open end of 
the test section and it is also the complement of x / d ( s / (T^'x/ d where ^ is the 
distance from the drive end of the test section). The Table 7.1 below gives the 
dimensionless operating parameters for a typical heater tube of the Stirling engine, the 
experimental test and the numerical simulation. 


Table 7.1 : Experimental and Numerical operating parameters. 


Parameters 

SPRE Heater 

Experiment Test 

Numerical Sim. 


11700 

11840 

11840 

Pa 

80.0 

80.2 

80.2 


1.03 

1.22 

1.22 

L/D h 

71.0 

60.0 

60.0 


1.2.2 Experimental Observations 

The preliminary experimental results carried out for different non dimensional 
parameters identified transition from laminar to turbulent flow by two mechanisms : 

► Convective triggering by the incoming turbulent fluid. 

► Instability of the developing boundary layer prior to the arrival of the convected 
turbulent fluid. 

Based on this observations a semi-empirical transition model has been proposed 
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(see Simon et. al.,1992) which accounts for the transition from laminar to turbulent flow 
effected by convective triggering of the incoming turbulent slug. The model is still in its 
seminal stages and hopes to improve the predictions for the skin friction and heat transfer 
coefficient by identifying the laminar and turbulent portions of the cycle at a given axial 
location. Further discussion of the model is beyond the scope of this study and can be 
found in Simon et al (1992). 

7.2.3 Numerical results and Comparison 

In order to correctly verify the numerical predictions care must be taken to 
accurately simulate the experimental conditions. Figures (7.11a) to (7.1 Id) shows how 
the friction factor compares against the experimental findings at the four axial locations, 
x / ^ 0.33,16,30,44 and the comparisons were made for the results obtained from the 
first oscillation cycle (i.e. the velocities have not settled to their quasi-steady states). The 
reason why the first oscillation cycle was chosen because the turbulent flow profiles, as 
was observed in the experiments, showed a cycle to cycle independence. That is, upon 
flow reversal the fluid was stagnant in the whole tube or had no history of the effects 
from previous cycle. A probable reason is due to the high mixing (turbulent diffusivity) 
capacities of turbulent flows, the velocity gradient near the wall were effectively 
smoothed out. Hence when the mean flow decelerates to zero velocity, the complete flow 
field in the tube also achieves no motion. Unlike in laminar flows, where due to relatively 
low diffusion the near wall velocity can be substantial compared to the core flow near 
flow reversal thus exhibiting a cyclic dependency. The experiments suggest a transient 
behavior in the sense that as the cycle begins slug (uniform inlet profile) fluid accelerates 


W 'll 

as a inviscid (flat velocity profile) flow close to the centerline of the tube having a 
growing boundary layer on the wall. These conditions are identical to when the fluid 
accelerates from rest in the beginning of first cycle even though the experimental data 
points are obtained after the flow has achieved a statistically steady state. 

During the laminar portion of the cycle the friction factors predicted numerically 
at the four axial positions compare very well with the experimental data as shown in the 
Figures (7.1 la)-(7.1 Id). The experimental data depart or show higher friction factors than 
the computed data after a particular velocity phase angle at each of the axial position 
except x / cF* 0.33. This phenomenon is attributed to the convective triggering of the 
boundary' layer from laminar to turbulent flow due to the arrival of the turbulent slug. 
Whereas at x/cTO. 33 which is close to the inlet the boundary layer is very thin and 
stable preserving its laminar state throughout the cycle, but as we go further down along 
the tube the growing boundary layer becomes very sensitive to the turbulent core and 
transitions to turbulent flow, resulting in higher skin friction coefficient. The reason for 
the high friction factor for the numerical results is that after the flow becomes turbulent 
the fully-developed turbulent flow correlations were used. This was done to test an 
empirical turbulence model (Simon et. al.,1992) for oscillatory flows to see how the 
predictions compare with the experimental result for the complete cycle. Also at 
x / D h ”0.33 the calculated friction factor are over predicted compared to the 
experimental data due to different inlet geometry' and being close to the inlet sensitive to 
the "entrance" effects. But as one goes down the tube "entrance" effects subside and 
predictions improve expectedly. 
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Figure 7.11 (a),(b),(c),(d). Comparison between Experimental and Numerical 
Moody friction factor plotted versus velocity phase angle at four different axial 
locations, namely, at x/D ■ 0.33,16.0,30.0,44.0 respectively. 
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Figure 7.11 (a),(b),(c),(d). Comparison between Experimental and Numerical 
Moody friction factor plotted versus velocity phase angle at four different axial 
locations, namely, at x/D m 0.33,16.030.0,44.0 respectively. 
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7.2.4 Brief Summary 

The experimental setup has been described in sufficient detail and a basis 
for comparison with the numerical computations derived. The rationale for using transient 
as opposed to quasi steady computed data was explained since the experimental data was 
collected after the flow inside the tube achieved a statistically steady state or no cyclic 
variations. Finally the accuracy of the code is corroborated by good agreement with the 
experimental data in the laminar portions of the cycle. 



7.3 Friction Factor and Entrance Effects 

The pressure drop required to set the fluid in motion is dependent on the so called 
friction factor ( c f ) which is a nondimensional wall shear stress (x w )- In the design of 
heat exchanger it is very important to optimize between the thermal and friction losses. 
The C{ ' s one suc h quantity which gives an estimate of the pumping power required to 

set the fluid in motion in an heat exchanger. The c f is defined as: 

c - (7.15) 

f (l/2)-p^i 2 

Here y . is the average instantaneous velocity across the cross section of the channel. 
Hence it shoots to infinity at the flow reversal instants of the cycle (O' 1 , 180 and 360 ) 
by the above definition since the average instantaneous velocity is zero at these points in 
a oscillation cycle. 

Figures 7.12, 7.13 and 7.14 plots the normalized friction factor versus the 
dimensionless axial distance at different velocity phase angles for the Cases C 2 , C 5 and 
H, respectively. Only curves for half oscillation cycle (0° to 180 p ) is shown (excluding 
the flow reversal points) due to the symmetry (in time and space) and the definition of 
r The friction factor is normalized by the instantaneous full developed steady state 
laminar friction factor correlation (.i.e. C f i?e i = 16 )• 1° l h' s particular half of the cycle 
the flow enters from the left end of the tube and exits at the right end of the tube ( x j D h 
- o to x / v h “ 60 or 70). From these figures the effect of Re mx ' Va and A r can be 
seen and these are described next. 

The and a together control the magnitude or level of the friction factor 


and any instant of the cycle. From Fig. 7.12 and 7.13 representing cases C 2 and C 5 one 
can see that the c value is higher for case C 5 than for case C 2 because of the Iow^ 
even though the pa is lower for case C^. This is can be attributed to the fluid inertia 
associated with the smaller fluid penetration into the tube, hence more pumping power 
has to be introduced in the case of low z - In contrast the heater case H 2 has lower 
values of friction factor because of the high ^ (see Fig. 7.14). Although for the same 
value of relative amplitude of fluid displacement (A r ) th e magnitude of the friction factor 
increases with p e . 

At some instants of the oscillation cycle the friction factor is negative, as can be 
seen fromFigures 7.12-7.14 at 150° velocity phase angle. Though it is unrealistic for the c f 
to be negative it simply means that the viscous forces augments the pressure forces. The 
negative value is caused by "backflow" at the walls or the velocities in a small viscous 
region close to the wall are flowing in a direction opposite to the primary flow direction. 
This situation is typical of high y a flows such as the cases C 2 , C 5 and H 2 , where the 
viscous effects are concentrated near a small region close to the wall (Stoke’s layer) with 
an inviscid core. The "backflow" or flow reversal at the wall occurs in the decelerating 
portion of the cycle and happens earlier for higher va which can be deduced from the 
Figures 7.12-7.14. It should be noted here that the "backflow" at the wall is characteristic 
of laminar oscillatory flows and does not exist when the flow is turbulent. The effect of 
flow reversal can also be noticed at the tube exit { x/ D h “ 60 or 70) during the 
accelerating portions of the cycle (30°) wherein the friction factor drops to value lower 
than the asymptotic value. Another effect of Valensi number ( y a ) is the amplitude of the 



friction factor at any axial location decreases with the y a (see Figures 7.12-7.14). The 
friction factor at any axial location is made up of a mean value plus a harmonic 
component, the amplitude here refers to the amplitude of the harmonic function. 

The "entrance" effects are directly related to the ^ value. Since it is difficult to 
define fully developed flow in an unsteady situation it is difficult to define an "entrance" 
length for oscillatory flows using the standard definition used for steady or unidirectional 
flows. But for low ^ such as case C 5 (Fig. 7.13) the friction factor is almost constant 
along the length of the tube at any given time except for a small region close to the 
entrance. And, this small length of the tube (the "entrance" length) where the friction 
factor drops from infinity to the constant asymptotic value, grows longer as the cycle 
advances in time. Therefore the "entrance" length behaves unsteadily and for the case C 5 
reaches a maximum value of 25 diameters ( x / D h " 25). Whereas for Cases C : and H ; 
(Figures 7.12 and 7.14) the friction factor does not reach to an asymptotic value at any 
instant of the cycle hence an "entrance" length is hard to define. The reason does not 
reach a asymptotic value can be ascribed to the high operating ^ values. But the drop 
in the value of friction factor from infinity to a lower value occurs over a small length 
for all the cases suggesting that the "entrance" effects are almost negligible. This 
observation is in contrast to steady flows where the "entrance" length increases with the 


Reynolds number. 




Figure 7.13. Normalized friction factor (C f ) versus the dimensionless axial 


distance at different velocity phase angles. 

Case C 5 :: tfe^-20000, Va-400 and A, -0.357. 
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Figure 7.14. Normalized friction factor ( C t ) versus the dimensionless axial 

distance at different velocity phase angles. 

Case H ; :: tfe^-16500, Ver 88 and A r -1.34. 
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CHAPTER VIII 

OSCILLATING FLOW HEAT TRANSFER 

Once the momentum and continuity equations are solved for then the energy 
equation can be solved for the thermal field. In this chapter the results of the thermal 
analysis for the individual components of the Stirling engine are presented. First, the 
thermal analysis for the regenerator which is a conjugate heat transfer type problem is 
presented. In the section the code validation aspect is also covered. Second, the results 
for the heater and cooler are presented with two different type of boundary conditions 
(Symmetric and Asymmetric Inflow) are discussed and presented. 

8.1 Conjugate Heat Transfer and Code Validation 

In this section the conjugate heat transfer phenomena occurring in the 
regenerator of the Stirling engine is discussed. A comparison of the numerical predictions 
with the analytical results for the temperature profile is also made. 
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8.1.1 Kurzweg Analysis (Analytical Solution) 

For the two-parallel-plates channel with solid walls and connected at the 
end to two reservoir at different temperatures (see Fig. 7.1), Kurzweg derived the solution 
for the temperature profile within the channel and the solid. This phenomenon is very 
similar to the one occurring in the regenerator of a Stirling engine. In the Stirling engine 
the regenerator is placed between the cooler and the heater which ensures a temperature 
gradient between the ends of the regenerator throughout the cycle. In the first half of the 
cycle the flow enters from the hot end and heat is absorbed by the regenerator and in the 
next half cycle when the flow enters from the cold end the absorbed heat is released to 
the cold fluid. Thus regenerator acts as a heat source or heat sink during a complete 
cycle. 

The temperature profile for along the channel for a fully developed velocity was 
derived by Kurzweg (1985a) for the geometry show in Figure 7.1. He assumed the 
temperature profile to be given by: 

T(x,y, t) = [T x + yag(ti) e ifi>t ] (8ll) 

Here the term r is the constant linear temperature gradient along the channel i.e. the 
temperature is assumed to be varying across the channel (normal j, direction) superim- 
posed on a constant axial temperature. And the function gr(Ti ) captures the variation 
in the normal direction and for the fluid region is given by g f (t| ) • 


g f { ti) = K x cosh ( fiPr a r) 


Xpe + — f ( t) ) 

a 4 Pr (Pr-l) a 2 (Pr-l) 


( 8 . 2 ) 


and for the solid portion is given by g s (r|) : 
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g s ( r| ) = K 2 cosh [ JToPra ( e - 1| ) ] 


The constants r and ^ are: 


-).Pe 

a 4 (Pr-1) Pr cosh WiPr a ] 


k y/Pr tanh [ Jig] tanh WicFra (€-1) ] | 
k tanh [v'TPra ] + /o tanh [JioPra (e-l) } J 


( 8 . 4 ) 


K x cosh [ s/i Pr a ] + 


*2 = 


A. Pe 

a* Pr ( Pr -l) 


cosh [ yj JoTr a ( e - 1 ) ] 


( 8 . 5 ) 


8.1.2 Numerical Predictions and Comparison 

Cases R, and R ; are the test cases used to compare the numerical predictions for 
the conjugate heat transfer problem with the Kurzweg analysis (1985a) presented above. 
The operating parameters for these cases are listed in Table 8.1. Since the solution for 
the temperature profile (Eqs. 8. 1-8. 5) are based on the flow being driven by a sinusoidally 
varying pressure gradient they have to be modified to account for the sinusoidal velocity 
input boundary conditions. 


Table 8.1 : Test cases investigated for the conjugate heat transfer. 


^9 


va 

L/D } 

Rfc stall 

T 

* west 

^east 










R, 

75 

2.5 

60 

n/a 

274 

273 

0.250 

R, 

12000 

400 

60 

n/a 

274 

273 

0.250 
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Essentially this involves two changes, first there will be a phase difference (<J> r ) 
added to the time variable ( u ,) and the function f ( n ) in Eq. (8.2) should be replaced 
by f (rj ) found from Eq. (7.8). The rationale behind this was explained in the fluid 
flow section. With these modifications the analytical temperature profile for a sinusoidal 
velocity input is given by: 


T{x,y, t) 


r t 

L x x 


yag(r\) 


^ ( i o> t + ) j 


(8.6) 


Where g (n ) is got from Eq. (8.2) or Eq. (8.3) depending on the fluid or the solid portion 
respectively. As mentioned earlier the function jf ( q ) should be replaced by the function ( t| ) 
given by eq. (7.8). The phase difference is given by: 


<|) r = 180° + 4> v • ’ 
The (b is gotten from eq. (7.6) and the 180° addition is due to the fact that the 
temperature gradient (y) used by Kurzweg is opposite in sign to the one used b> the 
numerical simulation. 

Figure 8.1 shows the temperature profile for case R„ Re^" 75 and Va = 7 - £ ’- 
The instantaneous temperature ( t ' T x ) or m( > re appropriately the instantaneous 
temperature fluctuation is plotted versus the normalized distance from the centerline to 
the wall ( y direction) at different velocity phase angles ( from 30° to 360° with 30" 
increment). The symbols are used for the analytical solution and the dotted line for the 
present work, the profiles were compared at the axial mid plane (£,/;) so that the 
entrance effects are negligible. Since the Valensi number ( va) is low the P r ° file exhibits 
the familiar parabolic shape as found in the steady state solution. Also the temperature 
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Figure 8.2. Comparison between analytical and numerical temperature 
profiles at different velocity phase angles for the conjugate heat transfer 
problem, Case R 2 12000 and Va-400. 

[Symbols : Analytical -Kurzweg ;; Dotted Lines : Numerical ] 
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gradient at the wall changes sign with the flow reversal or the wall heat flux is in phase 
with the mean or inlet velocity. The temperature fluctuation is almost zero in the solid 
because of the high heat capacity of the solid hence the constant % in Eq. (8.3) tends to 
zero with the result that (t| ) is zero for the whole cycle. Otherwise for low y a an( J 
low 0 (ratio of heat capacities of the solid and fluid) the temperature fluctuations within 
the solid is as high as in the fluid. 

Figure 8.2 shows a similar figure for case R 2 , J?e “ 12000 and va“ 400 which 
is a high y a case. Some interesting points to notice about this case are the presence of 
an inviscid core just like in the velocity plots due to the high frequency and the sharp 
temperature gradient at the wall which is out of phase with the incoming fluid. Also the 
temperature fluctuation is zero within the solid due to both the high heat capacity of the 
solid and the high y a . In fact one of the effects of high frequency on the temperature 
fluctuations in the solid (for a low heat capacity) is to bring the fluctuations down to zero. 
But due to both these factors present in this particular case it is difficult to isolate which 
has a more pronounced effect on the temperature distribution. Further the temperature 
fluctuation for case R 2 is more (*0.17) than for case R, (*0.035) which is a direct of the 
frequency or the presence of an inviscid core. 

Form Figures 8.1 and 8.2 one can clearly see that the agreement between the 
analytical solution and the numerical predictions are excellent. The reason for the poor 
agreement for the high y a case than the low y a can be attributed to the limitation of the 
machine in evaluating complex numbers algebra. Therefore special limiting procedures 
had to be performed on Eq. (8.2)-(8.5) so as to simplify them. 
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8.2 Symmetric Temperature Inflow 

The results for the heat transfer in circular tubes for a symmetric temperature 
inflow conditions with constant wall temperature are presented now. This type of heat 
transfer situation is close to what is found in the cooler and heater tubes of a Stirling 
engine. Table 8.2 lists all the cases investigated which cover a wide range of Re m&x -> Va 
and & . All the cases beginning with "C" stand for the cooler conditions and "H" for the 
heater conditions. In all the cases the flow enters the tube from the west end of the tube 
for the first half cycle (0° to 180'’), and then reverses and enters from the east end of the 
tube for the next half cycle (180° to 360°). In the case of symmetric inflow the flow 
enters the tube with the same constant temperature (enthalpy) for both half cycles. All 
the calculations were done for 52x52 grid and the runs were made on a Cray YMP 8/8128 
(sn 1040). Each run took approximately 1000s of CPU time for each cycle and at least 
3 cycles were needed for cyclic convergence or for the temperature profile to settle down 
from one cycle to the other. The computaional domain with the appropriate boundary 
conditions is sketched in Fig. (8.3), due to symmetry the computations were done only 
for half of the tube. 

8.2.1 Temperature Profiles 

Figures 8.4 and 8.5 are 3D plots for the nondimensional temperature profiles for 
the whole tube for two cooler cases with the same ^ . The nondimensional temperature 
( | ( T-T w ) / ( T in -T w ) | ) is plotted against the pipe radius and the axial distance at 
different velocity phase angles ( from 0° to 180° with 30° increment). Since the 
temperature inflow is symmetric the other half cycle is a mirror image of the events 
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Table 8.2\ Test cases investigated for symmetric temperature inflow. 


TEST 


va 

L/D, 

Twall 

T 

* west 

T 

x east 


CASE 








c, 

15000 

175 



340 

340 

0.714 

c 2 

30000 

350 

60 

325 

340 

340 

0.714 

C 3 

60000 

700 

60 

325 

340 

340 

0.714 

C 4 

10000 

200 

70 

325 

340 

340 

0.357 

c, 

20000 

400 

70 

325 


Ea 

0.357 

C 6 

40000 

800 

70 

325 

340 

340 

0.357 

H, 

8250 

44 

70 

650 

620 

620 

1.340 

h 2 

16500 

88 

70 

650 

620 

620 

1.340 

Hi 

33500 

176 

70 

650 

620 

620 

1.340 
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occurring in the first half. The direction of the flow is pointed by the thick bold arrow 
at any instant of time. 

Fig. 8.4 shows the cooler case C 2 ( J2e majf "' 30000, Va = 350, l/ D h = 60 and ^ = 
0.714 ), the flow enters with constant temperature (340 TC) and gets cooled due to the 
presence of a cold wall maintained at a constant temperature j* = 325 Tv. One can see 
as the cycle proceeds the thermal front advances into the tube but due to the high y a the 
profile across the radius does not develop into the familiar parabolic shape. The cooling 
process occurring is very complex due to the presence of cold and hot regions, at the 
beginning of the cycle (0°) the fluid is cooler at the entrance (for this part of the cycle) 
x/D, " 0 than at the exit x/D h ,s 60. The presence of these cold and hot "spots" 
effects the heat transfer mechanism in the tube as the flow proceeds along the tube. If 
one looks at the trace of the centerline ( j“0) throughout the half cycle (0 to 180 ), the 
temperature drops to a minimum at some axial location at any point in the cycle and then 
monotonically increases due to the presence of the hot spot ahead (in the axial direction). 
During the course of the cycle as the thermal front advances this temperature minimum 
moves further down the tube or closer to the exit and the being less than 1, this 
minimum temperature stays within the tube. Further due to the high y a the temperature 
gradient at the wall ( r - ± 0.05) is very steep suggesting a high wall heat flux. 

Fig. 8.5 shows a similar 3D plot for case C, ( 15000, y a - 175, j_, / “ 

60 and ^ = 0.714 ), the cooling mechanism is similar to the one described for case C 2 
except for few perceptible differences. The Valensi number (y a ) for this case being 
lower than case C 2 the hot temperature core is thinner than the previous case which makes 






Figure 8.3. Geometry for flow inside Circular tube showing the computaio- 
nal domain and Symmetric temperature inflow conditions. 
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Figure 8.4 Three-Dimensional normalized temperature plots at different velocity phase 
angles for the circular pipe symmetric temperature inflow problem: 

Case C 2 : Re > - 30000 and Va - 330. 
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the temperature profile to be closer to a parabola (e.g. see at 0° velocity phase angle at x j jj h 
=0). Also the temperature gradient of the wall is less steeper thus implying a lower wall 
heat flux. Figures 8.4 and 8.5 one can notice the effect of same a z by the maximum 
penetration distance of the thermal front into the tube which is the same, despite the wall 
heat flux being different. The temperature profiles are effected due to different y a . 

Figures 8.6 and 8.7 are 3D plots for cases H 2 ( 16500, y a = 88, l/ D h = 

70 and a = 1 .34 ) and C, ( Re - 20000, y a = 400, L/ D, ~ 70 and a m 0.357 ) 
respectively. Each case is for different a ar >d case H 2 is representative of the operative 
conditions of a heater in the Stirling engine where case C 5 was chosen to isolate the effect 
of the relative amplitude of fluid motion. Once again the nondimensional temperature is 
plotted against axial distance { x / D h ) a °d the radius of the tube ( r )■ In Fig. 8.6 which 
is the heater case H 2 , one can observe that as the cycle proceeds the thermal front has 
penetrated the whole axial distance and the cold "spot" which exists at the exit ( x / D h ~ 
70) is pushed out completely from the tube (at 90°). Whereas for the case C 5 shown in 
Fig. 8.7 the hot "spot" (since it is run as a cooler) exists throughout the half cycle. Also 
one can notice that lower the Valensi number ( y a ) the more parabolic the profile is ahead 
of the thermal front (or the inviscid hot or cold temperature core is thinner for low y a . 
8.2 .2 Contour Plots 

g-2. 

In this section the temperature contours for three of the cases listed in Table 
will be presented and discussed. These contour plots augments the discussion presented 
in the temperature profile section (see above). 

Figure 8.8 shows the temperature contours at different velocity phase angles for 
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case C 2 {a ~ 0.714). Again the contours are shown only for half of the cycle (O' to 1 80 ) 
due to symmetrical inflow conditions with the hot fluid entering the tube (at 7= 340 ”K) 
while the wall is maintained at a colder temperature (at 7= 325 °K). At the beginning 
of the cycle (0°) one can see the presence of hot fluid in the core of the tube (bigger at 
the east/left end of the tube) due to the history of hot fluid from the previous cycle. As 
the cycle proceeds during the acceleration portion (0° to 90") advances down the tube and 
continues to advance during the deceleration portion of the cycle (90° to 180") although 
at a slower rate. On the other hand the residual hot front from the previous cycle is 
retreating out of the tube during the acceleration portion of the cycle and disappears 
completely during the deceleration portion of the cycle (90° to 180"). The presence of 
two hot fronts (one close to the entrance and the other close to the exit), in most of the 
cycle, is attributed to the J?e max being less than 1. 

Figure 8.9 shows the temperature contours for the heater case H ; (^ r “ 1 -34). The 
contours are shown only for half a cycle (0 n to 180°) with cold fluid (heater) entering the 
tube ( t = 620 TC) while the wall being maintained at a hotter temperature ( 7- = 650 °K). 
In this case the cold front advances into the tube from the left end during the acceleration 
portion of the cycle and continues to advance during the deceleration portion of the cycle. 
Whereas the residual cold front present at the right end of the tube (at 0") is pushed out 
of the tube during the acceleration portion and disappears completely during the 
deceleration portion of the cycle. It is interesting to note that the core of cold fluid in the 
tube is thinner than the case C 2 (Fig. 8.8) due to lower va- 

Figure 8.10 shows the temperature contours for the cooler case C 5 (^ r “ 0.357). 


10 ? 


Again, due to symmetrica] inflow conditions the temperature contours are shown only for 
half a cycle. The hot fluid enters the tube at 7 -= 340 °K while the wall is maintained 
at a colder temperature (at 7 = 325 °K) so as to allow the incoming fluid to cool. Just 
like in the previous two cases the entering hot thermal front advances into the fluid during 
the acceleration portion and continues to advance but at a slower rate during the 
deceleration portion of the cycle. Further the residual hot front concentrated near the exit 
(right end) retreats out of the tube as the cycle advances but, contrary to the previous two 
cases, does not disappear completely from the tube at end of the half cycle. This can be 
attributed to the very low value of ^ (0.357). Despite the high yg the hot fluid core is 
thinner at the axial center of the tube than at the ends of the tube which is due to the very 
low Since the hot thermal front does not penetrate completely into the tube this 
means that around the axial center of the tube (away from the ends) the heat transfer 
mechanism is dominated by molecular conduction than by convection. This implies that 
heat transfer rate is very poor away from the ends. 


Figure 8.8. Temperature contours at different velocity phase angles for the symmetrical 
temperature inflow problem: Case Q Re„^y m 30000 and Va ” 350. 

[ A: 340.00, B: 339.99, C: 339.92, D:339.70, E: 339.10 

F: 337.50, G: 335.00, H: 330.00 , 1:327.50, J: 326.50 ] 


Figure 8.9. Temperature contours at different velocity phase angles for the symmetrical 
temperature inflow problem: Case H> Re rmx - 16500 and Va “ 88. 

[ A: 649.50, B: 640.00, C: 630.00, D: 625.00, E: 622.50 

F: 620.90, G: 620.30, H: 620.08 , I: 620.01, J: 620.00 ] 







8.2.3 Section Average Temperature (j> ) 

Figures 8.11-8.13 shows the normailzed section average temperature (r ) across 
the cross-section of the tube versus the dimensionless axial distance ( xfD h ) at different 
velocity phase angles for three cases. The section average temperature is a area weighted 
temperature across the tube and is defined as: 



Here ^a is an element of cross-section area, and a is the total cross-sectional area of 
the tube or channel. All the figures (8.11-8.13) are shown only for half of the oscillation 
cycle due to the symmetrical inflow. 

Figure 8.11 shows the normalizd x a versus xf D h for the cooler Case C, ( va 
=350) at different velocity phase angles. At the beginning of the cycle (0°) most of the 
fluid is hot at x f ]j h =60 and cold at x/ D h = 0 due to the hotter fluid leftover from the 
prvious half-cycle. As the flow accelerates the fluid cools rapidly due to the presence of 
the cold wall and the cold fluid present close to the entrance ( x/D h “0), and this 
temperature drop continues till it reaches a minimum value at certain axial locations (see 
curves for 30°, 60°, 90°) after which it increases monotonically down the tube. This 
minimum value keeps changing its axial location as the cycle goes on. But the level of 
temperature keeps dropping after the minimum temperature axial location and can be 
easily seen from the Figure 8.11 at x/ D h = 60. The displacement of the minimum y 
is dependent on value of a z too. From Figure 8.12 which is for case H ; ( va = 88) this 
minimum value moves out of the tube as the cycle advances due to a being greater than 
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Dimensionless Axial Distance (X/Dh) 


Figure 8.11. Normalized Section average temperature (r a ) versus the 

dimensionless axial distance ( x/D h ) at different velocity phase 
angles. Case C-, :: i?e m . v , '30(X)0, Va-350 and A *0.714. 



Dimensionless Axial Distance (X/Dh) 


Figure 8.12. Normalized Section average temperature (T a ) versus the 

dimensionless axial distance ( xf D h ) at different velocity phase 
angles. Case C 5 :: /fe^-20000, Ma-400 and A z -0.357. 
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Dimensionless Axial Distance (X/Dh) 

Figure 8.13. Normalized Section average temperature (r a ) versus the 

dimensionless axial distance ( xfD h ) at different velocity phase 
angles. 

Case H ; :: JJe^-16500, Ver 88 and A r -1.34. 
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one. The "minimum" temperature for the heater should actually be read as the maximum 
temperature as the normalized j’ is plotted instead of the true temperature. In the heater 
Case H-, the basic heat transfer mechanism is the same as described above for case C 2 but 
the direction of heat transfer reversed since the function is to heat the fluid. The Valensi 
number ( y a ) effects the amplitude or temperature fluctuation at any given axial location 
with time. The heat transfer mechanism behaves similarly for the same ^ . 

8.2.4 Bulk Temperature 

The use of the bulk temperature is very prevalent in the design of heat exchangers. 
Also known as the mixing cup temperature it is a velocity weighted temperature defined 
as: 


_ Jj>UTdA 

* P 

Here p m and u m are the mean density and velocity respectively. The mean velocity is 
given by the mathematical expression: 


U„ 



( 8 . 10 ) 


In oscillatory due to "backflow" or conuterflow at the walls at high y a this 
definition breaks down apart from the fact that the mean velocity goes to zero at 0°, 1 80" 
and 360° giving rise to unpysical values for the bulk temperature. Due to these 
discrepencies some researchers (Patankar and Oseid, 1992) have used the modulous of the 
velocity in Eq. (8.9) instead of the true velocity. With this change the Bulk temperature 
( t ) definition is now given as: 
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Dimensionless Axial Distance (X/Dh) 

Figure 8.14. Normalized Bulk temperature ( T b ) versus the dimensionless axial 


distance ( x/D h ) al different velocity phase angles. 
Case C 2 :: iJe^-SOOOO, Va“350 and A r -0.714. 



Dimensionless Axial Distonce (X/Dh) 

Figure 8.15. Normalized Bulk temperature ( T b ) versus the dimensionless axial 

distance ( x / D h ) at different velocity phase angles. 

Case C, :: ^e^-20000, Va“400 and A z -0.357. 
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Figure 8.J6. Normalized Bulk temperature ( T b ) versus the dimensionless axial 

distance ( x/ D h ) at different velocity phase angles. 

' Case H : :: Re^ x m 1 6500, Ver 88 and A f - 1.34. 
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Figures 8.14-8.16 plots the normalized bulk temperature profiles versus the 
dimensionless axial distance ( x/ D h ) at different velocity phase angles for the Cases C,, 
C 5 and H : respectively. Once again the Cases were chosen such that they represent 
different ^ values as the effects and transfer mechanisms are similar for the same^ r 
value. The behaviour of the y is similar to the trend of x a except for few differences. 
These differences can be observed by comparing Figures 8.13 and 8.16 which are for the 
same heater Case H ; , the amplitude of the bulk temperature ( j* ) is different from the 
amplitude of section average temperature ( j* ) fluctuation at any fixed axial location and 
even the magnitude of normalized value is different. This differences are of relatively 
minor consequence compared to the phase difference that arises between the and the 
wall heat flux ($J') and which will be shown later in the chapter. The presence of a 
phase difference is due to the definition of j* (Eq. 8.11) where the product of the local 
axial velocity ( jj) and temperature ( x ) is taken to evaluate x b • And, when the y a is 
high there is phase difference between the velocity and temperature such that their 
integrated product over a cross-sectional area is different, i.e a situation where the j* 
is zero can arise if jj and t are 90° out of phase. 



8.2.5 Wall heat Flux (a ") 
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The wall heat flux ( <%J') is the parameter which governs the effectiveness of ay 
heat exchanger. Figures 8.17-8.19 plots the q " normalized with the maximum inlet 
enthalpy ( qjj ct ) versus Wn at different velocity phase angles for the three 

Cases C 2 , C 5 and H 2 representing different ^’s. The negative wall heat flux values 
indicates that the direction of heat transfer is from the fluid to the wall, i.e the fluid gets 
cooled. The general trend is the wall heat flux increases along the length of the tube at 
instant of the cycle and then increases to a maximum or minimum value (depending upon 
whether it's a cooler or heater case) along the length of the tube and then starts to 
decrease or increase for the rest of the channel length. 

Since the wall heat flux is an important design parameter an attempt has been 
made to correlate the $jr as a function of A r , Re^ x , Va> and x/ D h for all the test 
cases investigated. The correlation process consisted of two important steps: 

► In the first the advantage of symmetric inflow conditions was taken. That is, the 
correlation effort was concentrated for only half oscillation cycle. 

► The second step involved Fourier analyzing the data for one half oscillation cycle 
(0° to 180°) and then correlate the resulting harmonics. 

The correlations were done using curve fitting software and was sufficiently tuned to 
handle multiple curve such as arising in oscillatory flows. The correlation arrived at has 
the form: 


A look at the above correlation reveals that a " consists of a mean value ( a" ) 

^ W \d]71 ' 



( 8 . 12 ) 
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= <r m ' f+ q"‘Sin{ 2 wt+<t> 2 ) 


Q it *"> ^2 

and an even harmonic component. <%" is the amplitude of the harmonic component and $ 2 

phase difference relative to the mean or inflow velocity and is a heat flux 

normalization factor used. As it can be seen the variation of wall heat flux is mathemati- 
cally and literally complex in oscilatory flows. All the coefficients are functions of all 
the nondimensional parameters and are mathematically found to be: 


fj \J V9 Q ^ ^ . -A a 

f! — v niax I -i c c 1 n -5 z . , 


i ' __ 


0 .312 -o . 1 -a; 


{ 1 + 5 . 5xl0' 5 -e ^(0 . 5~A r ) •( x/D h ) 


+ [0.84 + (1.05l//^)]-(A r -0.5Ml-A r )- v /3c7^ 


:8 . 13) 


+ [0.21+0.1437 *ln ( } 

The mean wall heat flux ( a" ) takes the following form: 

= (1.16) -[-1 +0.15 -(x/D h ) +7.5 7*7^] (8.14! 

and the other coefficients are: 

<J - 2 = -0.2+ 3 .0xicr 31 -e‘ <x/lv + [1/ (0 .26 + 0 .07 3A r ) ] • (x/D h ) ~ 0 - 5 (8.15) 

$ 2 = 262. 4 -24. 86-[ln(A r )] 2 -e [ ' 1 - e?4 ' 2 - 6A r 2 ' ln( ^ )1 -(x/D h ) 2 - 


( 1 + e [-4.75«-a.4-Vln(A r ,]. u/ ^ ) , (8.16) 

The above equation is valid for the first half oscillation cycle (0° to 180°) for the 
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Figure 8.17. Normalized Wall heat flux (&J 1 ) versus the dimensionless axial 
distance ( x / D h ) at different velocity phase angles. 

Case C 2 :: 30000, Ver 350 and A r ~ 0.714. 



Figure 8.18. Normalized Wall heat flux ( <%J ' ) versus the dimensionless axial 
distance ( x/D h ) at different velocity phase angles. 

Case C 5 :: 20000, Ver 400 and A x ~ 0.357. 
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Figure 8.19. Normalized Wall heat flux ( &J 1 ) versus the dimensionless axial 
distance ( x / D h ) at different velocity phase angles. 

Case H, :: Re^ m 16500, Ver 88 and A r - 1.34. 
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second half oscillation period the same equation is to be used but with x/D h being 
replaced by [l - x/ D h \ in all the above equations. 

It can be seen that the correlation coefficients take a very complicated functions 
of and x/ D h ■ The above equation was used to validate the numerically calculated 
wall heat fluxes <%J' • Figures 8.20-8.22 shows plots the ratio of numerical wall heat flux 

to the wall heat flux calculated using the correlation (Eq. 8.12) at select velocity phase 
angles for the cases C 2 , C 5 and H 2 . From the figures it can be concluded that the that the 
correlation predicts the numerical q " within 10% for the three cases for the velocity 
phase angles shown. The prediction is especially good near the centre of the pipe, the 
relatively poor prediction for low ^ values (see Fig. 8.23) has been caused by optimizing 
trying to optimize the correlation for higher ^ values (Cases C 2 and H 2 ). 

8.2.6 Heat Transfer Coefficient 

The heat transfer coefficient is usually represented by the Nusselt’s number ( jy u ) 
which is defined as: 


Nu 


[Ar] k 


( 8 . 17 ) 


Here the AT * s the reference temperature difference can be based either on section 
average temperature ( x a ) or the bulk temperature ( x b ) that is A T = [ T w ~ T a ] or [ j - 

Since the heat transfer coefficient is the variable used in most heat exchanger 
design, care should be taken as to how it is defined. As mentioned earlier the r and 
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Dimensionless Axial Distance (X/Dh) 

Figure 8.20. Ratio of Numerical Wall heat flux to the correlated wall heat flux 
versus the dimensionless axial distance ( x/ D h ) at different 
velocity phase angles. 

Case C 2 :: Re 30000, Ver 350 and 0.714. 
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Dimensionless Axiol Distance (X/Dh) 

Figure 8.21. Ratio of Numerical Wall heat flux to the correlated wall heat flux 
versus the dimensionless axial distance ( x/ D h ) at different 
velocity phase angles. 

Case C 5 :: Re^- 20000, Ver 400 and A r - 0.357. 






Figure 8.22. Ratio of Numerical Wall heat flux to the correlated wall heat flux 
versus the dimensionless axial distance ( x/ D h ) at different 
velocity phase angles. 

Case H 2 :: 16500, Va= 88 and A z “ 1.34. 



T b behave very different from each other in oscillatory flows thus making it very 
ambiguous as to how the temperature difference should be taken to evaluate the heat 
tansfer coefficient. This fact is demonstrated by the Figures 8.24a and 8.24b, in the 
former figure the Nusselt’s number based on x a while in the latter it is based on x b for 
the Case C 2 . When the section average temperature ( j* ) is used the the wall heat flux, 
Nusselt’s number and are all in phase with each other. Whereas when j* (Bulk 
temperature) is used they are all out of phase with each other (see Fig. 8.24b). Moreever, 
when x b is used the Nusselt’s number shoots to infinity at flow reversal points (0°, 180 n 
and 360°). It should be noted that the two plots show symmetry in time i.e. the thermal 
cycle is repeated twice for one flow cycle, this due to the inflow temperature symmetry 
and the plots are made at the channel mid-plane. 
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Figure 8.23a. Normalized values for Nusselt’s number ( Nu ), Wall 

heat flux ($„") and Temperature Difference (A T-T-T ) based on the 
section average temperature. 

Case C 2 :: Re^- 30000, Ver 350 and a z = 0.714. 



Figure 8.23b. Normalized values for Nusselt’s number ( Nu ), Wall 

heat flux (<$J') and Temperature Difference (A T- T b - T u ) based on the 
bulk temperature. 

Case C, :: Re^~ 30000, Ver 350 and A z ~ 0.714. 
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Figure 8.24a. Normalized values for Nusself $ number ( Nu ), Wall 

heat flux (<J J f ) and Temperature Difference (A 7 1 * r a - r w ) based on the 
section overage temperature. 

Case H 2 :: Re M - 16500, Va= 88 and A r - 1.34. 



Figure 8.24b. Normalized values for Nusselt’s number ( Nu ), Wall 

heat flux ( Qj 1 ) and Temperature Difference ( A T- T b - T w ) based on the 
bulk temperature. 

Case Hi :: Re max - 16500, Va- 88 and A z ~ 1.34. 




8.3 Asymmetric Temperature Inflow 

In this situation the flow enters the tube with a different temperature from one-half 
cycle to the other. In the Stirling engine heat exchangers this is the type of temperature 
boundary condition one encounters due to the fact the heat exchangers are arranged 
connected together. Table 8.3 list the cases investigated under the heater and cooler 
conditions as can be seen these are similar to the symmetric temperature inflow condition 
cases except for the inlet temperatures. 

The asymmetric temperature inflow heat transfer situation is complicated by the 
presence of an additional driving potential or the presence of axial heat conduction. 
Figure 8.25 shows an schematic representation of the cooler with asymmetric inflow 
temperature boundary conditions. In the first period (0° to 180") the hot fluid ( r ) 
enters from the left end of the tube and the second half period (180° to 360") it enters 
with a colder temperature ( T ) from the right end, with the wall being maintained at a 
colder temperature ( x ) than the either of the inflow temperature (i.e. T < T < 
T h )■ The heat transfer mechanism has two driving temperature potentials, ( y - x ) 
and ( T h - T c ) and depending upon which driving potential is greater the heat transfer 
mechanism is accordingly effected. In the present study, both ( r - r ) > ( r - 
X c ) and ( x c - T w ) < ( T h - T c ) have been investigated, from Table 8.3 it can be 
seen that the former condition exists for the heater Cases HA„ HA : and HA, and the 
latter for the cooler Cases CA,, CA-, and CA V The presence of axial driving potential 
plays a very important role especially in oscillatory flows where the axial heat transport 
is greatly augmented (See Kurzweg, 1985a) due to flow oscillation. And if the flow is 
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Figure 8.25. Asymmetrical temperature inflow boundary conditions 

representation for the heater. 
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at a high y a , additional complications can arise due to the backflow at the wall thus 
generating a additional temperature potential between the Stoke's layer and the central 
core of the fluid. 

Keeping the above observations in mind the temperature profiles for two 
representative Cases are discussed next. 

8.3.1 Temperature Profiles 

Figures 8.26 (a),(b) and 8.27 (a),(b) shows the 3D temperature profiles for the 
Cases CA 2 { Re ^ y =30000, y a = 350, =0.714, y^ = 350 % y^ = 330 Tv) 

and HA, (^ e *16500, y a = 88, A =1.340, y = 630 TC, y = 610 °K) at 

different velocity phase angles. In the discussion to be presented next the inflow' 
temperature from the west (left) end of the tube y will be referred as r (Hot 
end) and the inflow temperature from the east (right) end of the tube r will be 

u east 

referred to as x c (Cold end). Since the inflow is asymmetric with time the profiles are 
shown for the complete cycle. 

Figure 8.26 (a),(b) shows the 3D temperature profile plots for Case CA, ( y a = 
350) in the whole tube for every 30° increment of the oscillation cycle (30 o ,60°,90 o ,..., 
360°). In the first half cycle (0° to 180°) [ Fig. 8.26a ] hot fluid with temperature y = 
350 °K enters from the west end of the tube and in the next half cycle (180° to 360°) [ 
Fig. 8.26b ] the flow enters from the east end with a lower temperature y = 330 "K 
with the wall being maintained at the lowest temperature y = 325 °K. At the 
beginning of the cycle the flow enters from the hot end with x h (550 °K) with most of 


Figure 8.26a. Three-Dimensional Temperature plots at 
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the tube at a colder temperature due to the temperature history of the previous cycle. 
When the cycle begins most of the fluid in the tube is at a colder temperature (330 
°K) due to the temperature history effect of the previous cycle. Then the hot thermal 
front advances into the tube from the west end of the tube ( x/D h = 0) ss the cycle 
proceeds (60 n ) pushing the cold thermal front out of the tube ( X /D, =60), just like in 
the symmetric inflow situation. But due to the presence of the driving potential ( - 

T c ) the hot thermal front does not appear till late in the accelerating phase (60"). The 
effect of high y a can be seen by the presence of thick hot central core and due to ( t - 
T w ) < ( T h - T c ), the temperature of the thick central core never falls below T 
(330 °K) in the whole cycle despite the presence of a much colder wall ( r = 325 °K) 
Further the effect of this driving potential can be seen forward of the hot thermal front 
where the temperature of the core drops almost quasi-linearly to the temperature of the 
cold front (see Fig. 8.26a at 60°, 90°). As the hot thermal front advances into the tube in 
time the fluid ahead of it (near x/D h "60) gets cooled by the colder wall but due to the 
weaker driving potential ( x c ~ T w ) the temperature does not fall below T . When 
the flow reverses (after 180° velocity phase angle) and starts from the right end of the 
tube ( x / =60) most of the fluid in the tube is now at the hotter temperature r (350 
°K). Also the driving potential ( x h ~ T c ) now acts in a direction opposite to the flow 
direction, thereby the hot fluid present near the left end of the tube ( X /D, “0) cools 

rapidly as the second half cycle proceeds (210" to 360"). Hence it aids the cooling 
process done by the colder wall temperature. 

Figures 8.27a and 8.27b shows the 3D temperature profiles for Case HA2 


= 16500, Va = 88, Ar =1.340, = 630 °K, y^ = 610 °K and ^ = 650 °K. 

In this case the longitudinal temperature driving potential { T w ~ T ) is greater than the 
axial temperature driving potential ( T h ~ T c )• Hence the heat transport from the wall 
dominates over the axial heat transport. In the first half of the cycle (0° to 180°) the flow 
enters from the west or "hot" end with y^ = 630 °K and the next half cycle (180° to 
360°) the flow reverses and enters from the east or "cold" end. Figure 8.27a shows the 
temperature profiles for the first half cycle, as can be seen the effect of the advancing 
thermal front appears quickly at around 60° velocity phase angle. Prior to that the flow 
within the tube is cooler than the incoming fluid due to the temperature history of colder 
fluid entering from the "cold" end (see Fig. 8.27b at 360°). Because of high a the front 
starting from x/D h “0 advances into the tube and continues to advance as the cycle 
proceeds. Since the flow is in the same direction as the axial driving potential ( y - the 
whole tube gets heated up to more than the incoming fluid temperature at the end of the 
half cycle (180°). When the flow reverses (180° to 360°) thereby opposing this axial 
temperature driving potential the heating process is not as effective as for the previous 
half cycle (see Figure 8.27b) despite a high longitudinal temperature driving potential. 
Therefore at 360° one can observe that the flow actually "cools" at the exit (in this 
situation x/ D h “0) rather than increase in the temperature. 

Hence it can be concluded that the axial temperature driving potential plays a 
crucial role in the effectiveness of the heat exchanger during the whole complete cycle. 












134 


CHAPTER IX 
CONCLUSIONS 


A summary of the computational results for the flow and thermal analysis in three 
different components of the Stirling engine, namely. Regenerator, Cooler and Heater are 
presented next. The cases investigated are summarized in Table 6.1 and they represent 
the operating conditions of NASA Space Power Research Engine in terms of , Va - 
L J D h i ancl A z - 1° actual engine operating conditions, all the cases examined (except 
for Case R,) should go through the laminar/transition/turbulent flow regimes throughout 
the cycle. This study was focussed on the effects of oscillatory flow under laminar flow 
conditions with constant thermophysical properties. 

Cases R, and R 2 resemble the Regenerator and have been modeled using the conjugate 
heat transfer problem with a two-parallel -plates channel. Cases C,, C 2 and C 3 as well as 


H„ H 2 and H 3 are modeled using the circular pipe geometry with isothermal wall, 
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resembling the Cooler and Heater respectively. 

The conclusions from the computations are as follows: 

► The fluid flow and heat transfer for regenerator (Case R, ) are quasi-steady 
i.e. the velocity and temperature profiles are parabolic and in phase. For 
high y a (Case R 2 ) "backflow" or flow reversal takes place near the wall 
during the decelerating portions of the cycle resulting in flat temperature 
and velocity profiles in the core of the channel. In these situation the 
viscous effects are restricted to a small region next to the wall and is 
known as the Stoke’s layer. 

► The Cases R, and R ; were not only used to study the foil type regenerator 
but also to validate the numerical code by comparing the computational 
results with the analytical solution under similar operating conditions. The 
comparison was not trivial because the analytical work was based on 
oscillatory pressure boundary condition for an infinite channel, while the 
computational work was done for a finite channel with sinusoidally varying 
velocity boundary condition. The numerical results compared very well 
with the reformulated analytical results. 

► For oscillatory flows in parallel-plates-channel the wall shear stress ( x ) 
and pressure drop (Ap) are augmented by almost factors of 5 and 80 
respectively, as the va increases. Also they are out of phase with the 
velocity by 45° and 90° as the y a increases. 

► The instantaneous entrance length is controlled by the ^ value. For high;^ 



values (H,-H 3 & C,-C 3 ) it is difficult to define a fully developed flow 
conditions as the instantaneous friction factor does not settle with the axial 
distance. But in general the entrance effects are restricted to a small region 
into the channel. 

The heat transfer results show for the cooler and heater with symmetric 
temperature inflow conditions show that the heat transfer coefficient (at the 
channel mid-plane) goes through two cycles per each flow cycle. Also, at 
high va temperature profile is out of phase with the velocity profile. 
The heat transfer mechanism are controlled by the temperature "history" 
effects which reveals itself by the presence of "hot" and "cold” spots during 
the cycle. 

The usual definition of the heat transfer coefficient in the case of 
oscillatory flows is ambiguous and limited due to the way the temperature 
difference is determined. The common practice of using the bulk or 
mixing cup temperature ( T b ) as the reference for evaluating the 
temperature difference breaks down due to the flow reversal close to the 
wall during parts of the cycle for high ya> especially for laminar flows. 
This ambiguity was resolved in the study by using the absolute value of 
the velocity in the definition of bulk temperature and then evaluate the heat 
transfer coefficient. 

Using the section average temperature (y* ) as the reference temperature 
for the temperature difference, the wall heat flux, temperature difference. 
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and the Nusslet's number are all in phase with each other and out of phase 
with the velocity. The above quantities are all symmetric about J 80° at the 
channel axial mid-plane otherwise the shape is different from the first half 
cycle to the second half cycle. 

* Using the bulk temperature as the reference temperature the wall heat flux 
and the temperature difference are out of phase with each other 
(accordingly Nusselt’s number) and with the velocity. Also, the 
temperature difference passes through zero and accordingly, the Nusselt 
number shoots to infinity. 

* The asymmetric inflow temperature conditions has a drastic effect on the 
heat transfer mechanism. The two temperature driving potential (one axial 
and the other longitudinal) control the direction of the heat transfer 
mechanism. 

9.1 Scope for Further Research 

The primary objective of this research was to investigate the thermal field in a two 
dimensional simulation of the Stirling engine heat exchangers. The analysis presented in 
this study was done independently for each of the Stirling engine. The present study can 
be augmented by: 

Connecting the components together, i.e., the regenerator, cooler and heater in one 
direct two dimensional simulation to better understand the effects of oscillating 
flow in each component. 
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adding compressibility effects by introducing the two pistons found in the engine 
and thus facilitate a complete two dimensional engine simulation, 
introducing empirical turbulence model into the code to simulate the transition 
from laminar to turbulent flow and thus improve the predictions for engine losses. 
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APPENDIX B 


DERIVATION OF PRESSURE CORRECTION EQUATION 

In the step 3) of the SIMPLE algorithm described in section 5.5.1 it was 
mentioned that a pressure correction equation is needed to correct the velocities such that 
the continuity equation. This pressure correction equation is derived from the discretized 
continuity and momentum equation. The continuity equation (3.1) when discretized with 
respect to cartesian coordinate system yields: 


eq. (B.l) is nothing but a mass balance over a control volume, substituting for the mass 
fluxes for the respective CV faces from eq. (5.6), one gets: 

% + p e u e?>y ~ P w Uj>y + p n V n 6x - p s V s 6x =0 (B.2) 

Since the velocities U * and V* obtained from the solution of the momentum equations 
in steps 1 ) and 2) of the SIMPLE algorithm were based on a guessed initial pressure field, 
the continuity equation (B.2) will not be satisfied yielding a mass source S m : 

P 0 - P 

+ p e U e by - p„.C/j6y + p n V‘bx - p s V^bx - S w ( E -3) 

Here the indicates the newly found velocities from the momentum equations. To 
eliminate the mass source S m , velocity corrections need to performed to enforce mass 
corrections, then the velocity corrected continuity equation becomes by definition: 




pf~ Pp 
h t 


6 


v + p 9 <u; + U' a )by - p w (UC + U\.)by 


(E . 4 ) 


+ P n ( V n + V' n ) bx - p s (V s ‘ + F / s )5x = 0 
Subtracting eq. (B.3) from eq. (B.4) yields an equation for the velocity corrections, 

namely: 

p e u ' e t>y - P w u'j>y + p n v' n bx - 9 s v' s bx = -S m (B.5) 

Now the velocity corrections are related to the pressure corrections P' . The discretized 
momentum equations are used to couple the velocity and pressure. By linearizing the 
velocity corrections it can be related to the pressure corrections as: 

o'. - - rjo. <.P' E -P' p ) (E . 6 ) 

a p 

U 'w ~ ~ (-J-) * (P'p-P'w) (B . 7 ) 

a P 

V 'n " " (-p) n (P's-P'p) (B . 8 ) 

= " {P, P- P 's > (B . 9 ) 


Substituting eqs. (B.6)-(B.9) into equation (B.5) yields the so called pressure correction 
equation: 


a pP P ~ N + a pp' E + S + N ” S m (B.10) 

a complete description of the coefficients can be found in Perid and Scheuerer (1989). 

This pressure correction equation (B.10) has the same structure as the other 
discretized equations hence can be solved using the same matrix solver. The boundary 
conditions for the pressure correction equations are derived from the velocities at the 
boundary. The way this is achieved is by setting the velocity corrections at the 
appropriate CV boundary face to be zero. 

Since the pressure and velocity are nonlinearly coupled the SIMPLE algorithm 
diverges if no underrelaxation is employed, in the present code the pressure is corrected 
by the following equation: 

p;* = Pp + a. p P'p 


where a p is the underrelaxation factor for pressure and is usually in the range 0.2 -0.4. 



